slacky.eu

Codice: Seleziona tutto

/*
    The Resynthesizer - A GIMP plug-in for resynthesizing textures
    Copyright (C) 2000  Paul Francis Harrison
    Copyright (C) 2002  Laurent Despeyroux

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
   
   Original GIMP plug-in stripped for use in IDL by Adam Goodenough Feb, 2004
   
   */

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

#include <glib.h>

//Number of patches in leaf node of kd-tree
#define LEAF_SIZE 16

/* Radius of neighbourhood in pixels, originaly was 3 : arbitrarily set */
#define RADIUS_OF_NEIGHBOURHOOD_IN_PIXELS 3

/* Utility classes ===================================================== */

//Pixel channel type
typedef unsigned int Pixelel;

typedef int bool;
const bool true = 1, false = 0;

//Intersperse the bits from several bytes (used to sort colors
//  allowing a fairly primitive quantization)
/* Optimisation by Laurent Despeyroux */
uint entwine(Pixelel *data, int n_el) {
  int i,j,k,l;
  uint result = 0;
  for(i=0,k=1,l=0;i<8;i++,k*=2)
    for(j=0;j<n_el;j++,l++)
      if (data[j]&k)
        result += 1 << l;
  return result;
}

//Structure to hold an "entwine"ified pixels  
typedef struct {
  int x,y;
  uint twine;
} Twine_list_item;

int twine_list_item_sorter(const void *a, const void *b) {
  if ( ((Twine_list_item*)a)->twine < ((Twine_list_item*)b)->twine)
    return -1;
  else if ( ((Twine_list_item*)a)->twine > ((Twine_list_item*)b)->twine)
    return 1;
  else
    return 0;
}

typedef struct {
  int x,y;
} Coordinates;


/***************************/
/***************************/
/*** Entropy_accumulator ***/
/***************************/
/***************************/

//Structure+functions used to calculate the number of bits required to
//encode some normally distributed numbers
typedef struct {
  uint *x;
  uint *x_squared;
  int n_els, n;
} Entropy_accumulator;

Entropy_accumulator *entropy_accumulator_new(int n_els) {
  Entropy_accumulator *neo = g_new(Entropy_accumulator,1);
  neo->n_els     = n_els;
  neo->x         = g_new(uint,n_els);
  neo->x_squared = g_new(uint,n_els);

  return neo;
}

void entropy_accumulator_free(Entropy_accumulator *e) {
  g_free(e->x);
  g_free(e->x_squared);
  g_free(e);
}

void entropy_accumulator_calculate(Entropy_accumulator *e,
                                   double *bits,double *entropy) {
  int i;

  if (e->n <= 1) {
    *entropy = (double)(e->n_els)*sizeof(Pixelel)*8.0;
    *bits = e->n*(*entropy);
    return;
  }

  *entropy = *bits = 0.0;
    
  for(i=0;i<e->n_els;i++) {
    double mu    = (double)(e->x[i]) / e->n;
    double sigma, L;

    double sum_squares = e->x_squared[i]-e->x[i]*mu;

    if (sum_squares < e->n-1)
      sum_squares = e->n-1;
    
    sigma = sqrt( sum_squares / (e->n-1) ); 

    L =
      e->n/2.0*log(2*3.1415926535) + e->n/2.0*log(sigma*sigma) +
      1.0/(2.0*sigma*sigma)*sum_squares;

    *bits += log(256)+log(sigma)+L+
             0.5*log((double)(e->n)*e->n/(2.0*pow(sigma,6.0)))
             -0.5*log(12)+0.5;

    *entropy += 0.5*log(2*3.14159265)+0.5+log(sigma);
  }

  //Convert from nits to bits
  *bits /= log(2);
  *entropy /= log(2);
}

void entropy_accumulator_clear(Entropy_accumulator *e) {
  int i;
  for(i=0;i<e->n_els;i++) {
    e->x[i]         = 0;
    e->x_squared[i] = 0;
  }
  e->n = 0;
}

void entropy_accumulator_add(Entropy_accumulator *e, Pixelel *data) {
  int i;
  for(i=0;i<e->n_els;i++) {
    e->x[i]         += data[i];
    e->x_squared[i] += (uint)(data[i])*(uint)(data[i]);
  }
  e->n++;
}

/**************/
/**************/
/*** Bitmap ***/
/**************/
/**************/

//Bitmap class with three dimensions (width, height, number of channels)
typedef struct {
  int width, height, depth, element_size;
  void *data;
} Bitmap;

Bitmap *bitmap_new(int element_size) { 
  Bitmap *neo = g_new(Bitmap,1);
  neo->element_size = element_size;
  neo->data = 0;

  return neo;
}

void bitmap_free(Bitmap *b) {
  g_free(b->data);
  g_free(b);
}

void bitmap_size(Bitmap *b, int w,int h,int d) {
  b->width = w; 
  b->height = h; 
  b->depth = d;
    
  g_free(b->data);
  b->data = g_malloc(w*h*d*b->element_size);
}

#define bitmap_at(b,type,x,y) \
  ( (type*)((b)->data)+((x)+(y)*(b)->width)*(b)->depth )


//Structure used to hold information about pixels in the output
typedef struct {
  bool use;
  bool done;
  int priority;

  int heap_index;
} Status;

/************/
/************/
/*** Heap ***/
/************/
/************/

//Heap used to prioritize the selection of pixels in the output
typedef struct {
  Status **data;
  int n_els;
} Heap;

void heap_put(Heap *h, Status *item,int index) {
  h->data[index] = item;
  if (item)
    item->heap_index = index;
}

Heap *heap_new(Status *items,int n_items) {
  int i;
  Heap *h = g_new(Heap,1);
  h->n_els = n_items;
  h->data  = g_new(Status*,n_items);
  
  for(i=0;i<n_items;i++) 
    heap_put(h,items+i,i);

  return h;
}

void heap_free(Heap *h) {
  g_free(h->data);
  g_free(h);
}

int heap_parent(int i) { return (i+1)/2-1; }
int heap_left_child(int i) { return (i+1)*2-1; }
int heap_right_child(int i) { return (i+1)*2; }

void heap_promote(Heap *h, Status *item) {
  int i = item->heap_index;
  int heap_parent_i ;
  while(i > 0 && 
          (!h->data[heap_parent_i=((i+1)/2-1)] || 
      h->data[heap_parent_i]->priority < item->priority))
    {
  /*  heap_put(h, h->data[heap_parent_i],i); */   /* the 2 next lines are the transcription of this call inline L.D. */
    if ((h->data[i] = h->data[heap_parent_i]))
      h->data[heap_parent_i]->heap_index = i;

    i = heap_parent_i;
    }
  heap_put(h, item, i);
}

void heap_become_heapish(Heap *h) {
  int i;
  for(i=0;i<h->n_els;i++)
    if (h->data[i]->done)
      h->data[i] = 0;
    else
      heap_promote(h, h->data[i]);
}

Status *heap_pop_top(Heap *h) {
  Status *top = h->data[0];
  int heap_right_child_i ;
    
  int i=0;
  for(;;) {
    int next = heap_left_child(i);
    if (next >= h->n_els) {
      /*heap_put(h, 0,i);*/   /* the next line is the transcription of this call inline L.D. */
      h->data[i] = 0;
      break;
    }

    if ((heap_right_child_i=heap_right_child(i)) < h->n_els &&
        h->data[heap_right_child_i] != NULL &&
        (!h->data[next] || 
    h->data[next]->priority < h->data[heap_right_child_i]->priority))
      next = heap_right_child_i;
      
   /* heap_put(h,h->data[next],i); */   /* the 2 next lines are the transcription of this call inline L.D. */
    if ((h->data[i] = h->data[next]))
      h->data[next]->heap_index = i;
    i = next;
      
    if (!h->data[i])
      break;
  }

  return top;
}

/******************/
/******************/
/*** KDTREE     ***/
/******************/
/******************/

// k-d trees are a way to do fast(ish) closest vector
// searches

typedef struct {
  int offset_no;
  Pixelel value;
} Pattern;

typedef struct {
  int x,y,el, offset;
  uint weight;

  //So we can sort from largest to smallest weight 
  // (provides a good speedup as distance calculations can often be aborted earlier)
  //bool operator < (const Offset &other) const { return weight > other.weight; }
} Offset;

int offset_sorter(const void *a, const void *b) {
  if ( ((Offset*)a)->weight > ((Offset*)b)->weight)
    return -1;
  else if ( ((Offset*)a)->weight < ((Offset*)b)->weight)
    return 1;
  else
    return 0;
}

/************/
/************/
/*** Tree ***/
/************/
/************/

struct _Node;

typedef struct {
  struct _Node *node;
  int best_dist;

  //bool operator < (const Queue_item &other) const { 
  //  return best_dist > other.best_dist; 
  //}
} Queue_item;


// Globals to avoid passing around these things incessantly 

// Array of Offset
GArray *Offsets;

// Array of Pattern
GArray *Patterns;

//More globals... array of Queue_item
GArray *node_queue;
int node_queue_count;

int cutoff_dist;
Pixelel *best;

void set_up_tree_globals() {
  Offsets    = g_array_new(FALSE,FALSE,sizeof(Offset));
  Patterns   = g_array_new(FALSE,FALSE,sizeof(Pattern));
  node_queue = g_array_new(FALSE,FALSE,sizeof(Queue_item)); 
}

void clean_up_tree_globals() {
  g_array_free(Offsets,TRUE);
  g_array_free(Patterns,TRUE);
  g_array_free(node_queue,TRUE);
}

void node_queue_upheap(int i) {
  Queue_item item = g_array_index(node_queue,Queue_item,i);

  for(;;) {
    int parent;
    
    if (i == 0)
      break;
      
    parent = heap_parent(i);
    if (g_array_index(node_queue,Queue_item,parent).best_dist <
        item.best_dist)
      break;

    g_array_index(node_queue,Queue_item,i) =
      g_array_index(node_queue,Queue_item,parent);

    i = parent;
  }
    
  g_array_index(node_queue,Queue_item,i) = item;
}

void node_queue_enqueue(Queue_item item) {
  if (node_queue_count == node_queue->len)
    g_array_set_size(node_queue,(node_queue->len?node_queue->len*2:1));

  g_array_index(node_queue,Queue_item,node_queue_count++) = item;
  node_queue_upheap(node_queue_count-1);
}

void node_queue_pop_top() {
  int i = 0;
  for(;;) {
    int child,right_child;

    child = heap_left_child(i);
    if (child >= node_queue_count)
      break;

    right_child = heap_right_child(i);
    if (right_child < node_queue_count &&
        g_array_index(node_queue,Queue_item,right_child).best_dist <
        g_array_index(node_queue,Queue_item,child).best_dist)
      child = right_child;

    g_array_index(node_queue,Queue_item,i) =
      g_array_index(node_queue,Queue_item,child);

    i = child;
  }

  g_array_index(node_queue,Queue_item,i) =
    g_array_index(node_queue,Queue_item,--node_queue_count);

  node_queue_upheap(i);
}

typedef struct {
  Pixelel min, max;
} Bound;

struct _Node {
  gboolean is_leaf;
  Bound *bounds;

  union {
    struct {
      Pixelel **patches;
      int n_patches;
    } leaf;
    struct {
      struct _Node *left_child, *right_child;
    } inner;
  } data;
};

typedef struct _Node Node;

void node_free(Node *node) {
  g_free(node->bounds);

  if (node->is_leaf)
    g_free(node->data.leaf.patches);
  else {
    node_free(node->data.inner.left_child);
    node_free(node->data.inner.right_child);
  }

  g_free(node);
}

void node_add_to_queue(Node *n)
  {
  Queue_item item;
  int i, dist = 0;
  for(i=0;i < Patterns->len;i++)
    {
    int offset_no = g_array_index(Patterns,Pattern,i).offset_no;
      
    if (g_array_index(Patterns,Pattern,i).value < n->bounds[offset_no].min)
      {
      dist += g_array_index(Offsets,Offset,offset_no).weight * 
        (n->bounds[offset_no].min-g_array_index(Patterns,Pattern,i).value);
      if (dist >= cutoff_dist)
   return;
      }
      else
        if (g_array_index(Patterns,Pattern,i).value > 
               n->bounds[offset_no].max)
          {
          dist += g_array_index(Offsets,Offset,offset_no).weight * 
            (g_array_index(Patterns,Pattern,i).value-n->bounds[offset_no].max);
          if (dist >= cutoff_dist)
            return;
          }
    } /* for i < Patterns->len */

  item.node = n;
  item.best_dist = dist;
  node_queue_enqueue(item);
  }

void node_do_it(Node *n) 
  {
  if (n->is_leaf)
    {
    int i,j;
  
    for(i=0;i<n->data.leaf.n_patches;i++)
      {
      int dist = 0;
      for(j=0;j<Patterns->len;j++)
        {
   int offset_no = g_array_index(Patterns,Pattern,j).offset_no;
   
        dist += g_array_index(Offsets,Offset,offset_no).weight * 
     abs((int)(n->data.leaf.patches
             [i][g_array_index(Offsets,Offset,offset_no).offset]) - 
         (int)(g_array_index(Patterns,Pattern,j).value));

   if (dist > cutoff_dist)
     goto skip;
        }

      best = n->data.leaf.patches[i];
      cutoff_dist = dist;

      skip: ;
      }
  
    }
    else
    {
    
    node_add_to_queue(n->data.inner.left_child);
    node_add_to_queue(n->data.inner.right_child);
  
    }
  }

Pixelel *query_tree(Node *tree) {
  //Clear the pq
  node_queue_count = 0;

  cutoff_dist = 1<<30;

  node_do_it(tree);

  while(node_queue_count) {
    Queue_item item = g_array_index(node_queue,Queue_item,0);

    if (item.best_dist > cutoff_dist)
      break;

    node_queue_pop_top();
    node_do_it(item.node);
  }

  return best;
}

Node *build_tree(Pixelel **items,int n_items) {
  int i,j;
  Node *node = g_new(Node,1);
  
  node->is_leaf = n_items <= LEAF_SIZE; 

  if (node->is_leaf) {
    node->data.leaf.patches = g_new(Pixelel*,n_items);
    node->data.leaf.n_patches = n_items;
    for(i=0;i<n_items;i++)
      node->data.leaf.patches[i] = items[i];
  }

  node->bounds = g_new(Bound,Offsets->len);

  for(i=0;i<Offsets->len;i++) {
    node->bounds[i].min = 255;
    node->bounds[i].max = 0;
  }

  for(i=0;i<n_items;i++)
    for(j=0;j<Offsets->len;j++) {
      int offset = g_array_index(Offsets,Offset,j).offset;

      if (items[i][offset] < node->bounds[j].min)
        node->bounds[j].min = items[i][offset];
      if (items[i][offset] > node->bounds[j].max)
        node->bounds[j].max = items[i][offset];
    }

  if (!node->is_leaf) {
    int best = 0;
    int best_offset, lower, upper, start_upper, pivot;
    Pixelel *temp;

    for(i=1;i<Offsets->len;i++)
      if (g_array_index(Offsets,Offset,i).weight*
            (node->bounds[i].max-node->bounds[i].min) >
          g_array_index(Offsets,Offset,best).weight*
       (node->bounds[best].max-node->bounds[best].min))
   best = i;

    best_offset = g_array_index(Offsets,Offset,best).offset;

    lower=0; 
    upper=n_items-1;
    for(;;) {
      Pixelel p_val = items[lower][best_offset];
      start_upper = lower+1;

      for(i=lower+1;i<=upper;i++)
        if (items[i][best_offset] < p_val) {
     Pixelel *temp = items[i];
     items[i] = items[start_upper];
     items[start_upper] = temp;
     start_upper++;
   }

      pivot = start_upper-1;
      temp = items[lower];
      items[lower] = items[pivot];
      items[pivot] = temp;

      if (pivot == n_items/2)
        break;

      if (pivot < n_items/2)
        lower = pivot+1;
      else
        upper = pivot-1;
    }

    node->data.inner.left_child = 
      build_tree(items,n_items/2);
    node->data.inner.right_child =
      build_tree(items+n_items/2,n_items-n_items/2);
  }

  return node;
}


/******************/
/******************/
/*** Parameters ***/
/******************/
/******************/

typedef struct {
  bool h_tile, v_tile;
  
  gint32 corpus_id, input_map_id, output_map_id;

  double random_weight, map_weight;
} Parameters;

bool wrap_or_clip(Parameters *parameters, Bitmap *image, int *xx,int *yy)
{ /* optimisation by Laurent Despeyroux, not important */
  if (*xx < 0)
    {
    if (parameters->h_tile)
      *xx += image->width;
    else
      return false;
    }
  else
    {  
    if (*xx >= image->width) 
      {
      if (parameters->h_tile)
        *xx -= image->width;
      else
        return false;
      }
    }
  
  if (*yy < 0)
    {
    if (parameters->v_tile)
      *yy += image->height;
    else
      return false;
    }
  else
    {
    if (*yy >= image->height)
      {
      if (parameters->v_tile)
        *yy -= image->height;
      else
        return false;
      }
    }

  return true;
}


    /*  ------------------------------------------------------------------  */
    /*                                                                      */
    /*  Here is the most important function : the one that does the work !  */
    /*                                                                      */
    /*  ------------------------------------------------------------------  */

static void run(void)
{ 

  int i,j,k,l,x,y;
  Parameters parameters;
  int input_bytes, map_bytes, map_pos, random_bytes, random_pos, bytes,
      n_pixels, n_corpus_pixels, patch_radius;    
  Bitmap *status, *weights;
  GArray *neighbours;
  Heap *heap;
  int weight_sum_len;
  double *weight_sum;
  Node *tree;

  /* ----begin---- */
  srand(time(0));
  
  /* Setup parameters */
  parameters.corpus_id = 2;
  parameters.input_map_id = -1;
  parameters.output_map_id = -1;
  parameters.v_tile = false;
  parameters.h_tile = false;
  parameters.random_weight = 0.25;
  parameters.map_weight    = 0.75;

  input_bytes  = 3;
  map_bytes   = 0;
  map_pos     = input_bytes;
  random_bytes = 1;
  random_pos   = map_pos + map_bytes;
  bytes        = random_pos + random_bytes;

  status = bitmap_new(sizeof(Status));          
                        //Have done this output pixel yet?, pixel priority
  weights = bitmap_new(sizeof(int));         
                        //Weighting given to neighbours when a pixel is added 

  patch_radius = RADIUS_OF_NEIGHBOURHOOD_IN_PIXELS ;
  
  neighbours = g_array_new(FALSE,FALSE,sizeof(Coordinates));

  /* Choose which neighbours to use */

  { 
    int radius;
    for(radius=1;radius<=patch_radius;radius++) {
      for(i=0;i<radius*2;i++) {
        Coordinates c;
        c.x = -radius+i;
        c.y = radius;
        g_array_append_val(neighbours,c);
        c.x = radius;
        c.y = radius-i;
        g_array_append_val(neighbours,c);
        c.x = radius-i;
        c.y = -radius;
        g_array_append_val(neighbours,c);
        c.x = -radius;
        c.y = -radius+i;
        g_array_append_val(neighbours,c);
      }
    }
  }
  
  /************MYMASK************/
  
  Bitmap *mymask = bitmap_new(sizeof(Pixelel));
   bitmap_size(mymask, 32, 32, 1);

   int xi,yi;
   for(yi=0;yi<mymask->height;yi++)
   for(xi=0;xi<mymask->width;xi++) {
      if ( xi >= 12 &&
         xi <= 23 &&
         yi >= 12 &&
         yi <= 19){
      ((Pixelel*)mymask->data)[yi*mymask->width+xi] = 255;
      }
      else ((Pixelel*)mymask->data)[yi*mymask->width+xi] = 0;
   }
   
    /************MYCORPUS************/
  
  Bitmap *mycorpus = bitmap_new(sizeof(Pixelel));
   bitmap_size(mycorpus, 32, 32, 4);

   for(yi=0;yi<mycorpus->height;yi++)
   for(xi=0;xi<mycorpus->width;xi++) {
      if ( xi >= 12 &&
         xi <= 23 &&
         yi >= 12 &&
         yi <= 19){
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4] = 255;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+1] = 128;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+2] = 64;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+3] = 0;
      }
      else {
      if (xi % 8 < 4){
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4] = 255;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+1] = 128;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+2] = 64;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+3] = 0;
      } else {
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4] = 64;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+1] = 128;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+2] = 255;
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+3] = 0;
      }}

   }
   
   /************MYDATA************/
   
   Bitmap *mydata = bitmap_new(sizeof(Pixelel));
   bitmap_size(mydata, 32, 32, 4);
   for(yi=0;yi<mydata->height;yi++)
   for(xi=0;xi<mydata->width;xi++) {
      ((Pixelel*)mydata->data)[(yi*mycorpus->width+xi)*4] = 
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4];
      ((Pixelel*)mydata->data)[(yi*mycorpus->width+xi)*4+1] = 
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+1];
      ((Pixelel*)mydata->data)[(yi*mycorpus->width+xi)*4+2] = 
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+2];
      ((Pixelel*)mydata->data)[(yi*mycorpus->width+xi)*4+3] = 
         ((Pixelel*)mycorpus->data)[(yi*mycorpus->width+xi)*4+3];
   }
  

  /* Allocate memory for output and status tracking */

  {

    bitmap_size(status,mydata->width,mydata->height,1);
    heap = heap_new(status->data,status->width*status->height);

    for(i=0;i<status->width*status->height;i++) {
      ((Status*)status->data)[i].use  = true;
      ((Status*)status->data)[i].done = (((Pixelel*)mymask->data)[i] != 255);
      ((Status*)status->data)[i].priority = 0;
    }

  }
    
  n_pixels = 0;
 
  {
    Pixelel *ptr = mydata->data;
    for(i=0;i<mydata->width*mydata->height;i++,ptr+=mydata->depth) {
      if (!((Status*)status->data)[i].done) {
        n_pixels++;
        for(j=0;j<input_bytes;j++)
          ptr[j] = 128;
      }
      ptr[random_pos] = rand();
    }
  }

  {
    Pixelel *ptr = ((Pixelel*)mycorpus->data) + random_pos;
    for(i=0;i<mycorpus->width*mycorpus->height;i++,ptr+=mycorpus->depth)
      *ptr = rand();
  }

  /* Work out which pixels from the corpus to use */

  {
    Pixelel *ptr;
    ptr = mymask->data;
    for(i=0;i<mymask->width*mymask->height;i++,ptr++)
      *ptr = (*ptr ? 0 : 1);
  }
  
  for(y=patch_radius;y<mycorpus->height-patch_radius;y++)
    for(x=patch_radius;x<mycorpus->width-patch_radius;x++)
      if (bitmap_at(mymask,Pixelel,x,y)[0])
        for(i=0;i<neighbours->len;i++) {
     Coordinates *n = &g_array_index(neighbours,Coordinates,i);
          if (bitmap_at(mymask,Pixelel,x+n->x,y+n->y)[0])
            bitmap_at(mymask,Pixelel,x+n->x,y+n->y)[0]++;
   }

  n_corpus_pixels = 0;
  
  {
    Pixelel *ptr;
    ptr = mymask->data;
    for(i=0;i<mymask->width*mymask->height;i++,ptr++)
      if (*ptr == neighbours->len+1) {
        *ptr = (*ptr == neighbours->len+1 ? 1 : 0);
        n_corpus_pixels++;
      } else
        *ptr = 0;
  }

  /* Sanity check */

  if (!n_corpus_pixels || !n_pixels) {
    if (!n_corpus_pixels)
      printf("The input texture is too small.");
    else
      printf("The output image is too small.");
   exit(1);
  }


  /* Examine neighbouring pixel interactions */

  bitmap_size(weights, mycorpus->width,mycorpus->height,neighbours->len);

  weight_sum_len = neighbours->len;
  weight_sum = g_new(double,weight_sum_len);
  for(i=0;i<weight_sum_len;i++)
    weight_sum[i] = 0.0;
    

  {
    double total_bits, total_entropy;
    Twine_list_item *sorted_list;

    uint mask;
    int start;

    Entropy_accumulator *accumulator = entropy_accumulator_new(input_bytes);

    entropy_accumulator_clear(accumulator);
    for(i=0;i<mycorpus->width*mycorpus->height;i++)
      if (((Pixelel*)(mymask->data))[i])
        entropy_accumulator_add(accumulator,(Pixelel*)mycorpus->data+i*mycorpus->depth);

    entropy_accumulator_calculate(accumulator, &total_bits, &total_entropy);
    
    sorted_list = g_new(Twine_list_item,n_corpus_pixels);
    
    for(i=0,y=0;y<mycorpus->height;y++)
      for(x=0;x<mycorpus->width;x++) 
        if (bitmap_at(mymask,Pixelel,x,y)[0]) {
     sorted_list[i].x = x;
     sorted_list[i].y = y;
     sorted_list[i].twine = entwine(bitmap_at(mycorpus,Pixelel,x,y), input_bytes);
     i++;
        }

    qsort(sorted_list,n_corpus_pixels,sizeof(Twine_list_item),
     twine_list_item_sorter);

    for(i=0;i<neighbours->len;i++) {
      /* Find best quantization level */

      int best_level = 0;
      double best_bits = 1e30;

      for(j=0;j<=(int)(sizeof(Pixelel)*8);j++) {
   double this_bits;

        mask = (uint)(-1);
   mask <<= j*input_bytes;
   
   entropy_accumulator_clear(accumulator);
        start = 0;
   this_bits = 0.0;
        for(k=0;;k++) { 
     if (k == n_corpus_pixels ||
         (sorted_list[start].twine&mask) != (sorted_list[k].twine&mask)) {
       double bits, entropy;
            entropy_accumulator_calculate(accumulator, &bits, &entropy);
       this_bits += bits;
            entropy_accumulator_clear(accumulator);
       start = k;
     }

     if (k == n_corpus_pixels)
       break;
       
     entropy_accumulator_add(accumulator,
       bitmap_at(mycorpus,Pixelel,
         sorted_list[k].x+g_array_index(neighbours,Coordinates,i).x,
         sorted_list[k].y+g_array_index(neighbours,Coordinates,i).y));
   }

        //Kinda dodgy, its possible the length won't be a... whats the
   // technical term... one hump function
   if (this_bits > best_bits)
     break;

   best_level = j;
   best_bits = this_bits;
      }
      
     /* Calculate weights */

      mask = (uint)(-1);
      mask <<= best_level*input_bytes;
   
      entropy_accumulator_clear(accumulator);
      start = 0;
      for(k=0;;k++) { 
        if (k == n_corpus_pixels ||
       (sorted_list[start].twine&mask) != (sorted_list[k].twine&mask)) {
     int weight;
     double bits, entropy;
          entropy_accumulator_calculate(accumulator, &bits, &entropy);

     weight = (int)( (total_entropy - entropy) * 65536.0 );
     if (weight < 0)
       weight = 0;
     
     weight_sum[i] += weight * (k-start);

          for(l=start;l<k;l++)
       bitmap_at(weights,int,sorted_list[l].x,sorted_list[l].y)[i] = 
         weight;
     
     entropy_accumulator_clear(accumulator);
     start = k;
   }
   
   if (k == n_corpus_pixels)
     break;
       
   entropy_accumulator_add(accumulator,
     bitmap_at(mycorpus,Pixelel,
       sorted_list[k].x+g_array_index(neighbours,Coordinates,i).x,
       sorted_list[k].y+g_array_index(neighbours,Coordinates,i).y));
      }
      
    }

    g_free(sorted_list);
    entropy_accumulator_free(accumulator);
  } /* Calculating weights */



  set_up_tree_globals();

  {
    double sum, map_sum;
    Offset offset;

    for(y=0;y<mycorpus->height;y++)
      for(x=0;x<mycorpus->width;x++) 
        if (bitmap_at(mymask,Pixelel,x,y)[0]) {
          int *weight = bitmap_at(weights,int,x,y);

     int sum = 0;
     for(i=0;i<neighbours->len;i++) {
       sum += weight[i];
     }

     if (sum == 0)
       for(i=0;i<neighbours->len;i++)
         weight[i] = 65536 / neighbours->len;
     else
       for(i=0;i<neighbours->len;i++)
         weight[i] = (int)(weight[i] * 65536.0 / sum);
        }

    for(i=0;i<weight_sum_len;i++) {
      weight_sum[i] = pow(2.0,weight_sum[i]
        /n_corpus_pixels
      /65536.0);
    }
    
    sum = 0.0;
    for(i=0;i<neighbours->len;i++)
      sum += weight_sum[i];
    
    map_sum = 0.0;

    for(i=0;i<neighbours->len;i++) {
      offset.weight = (int)(weight_sum[i]*65536.0 / sum);
      offset.x  = -g_array_index(neighbours,Coordinates,i).x;
      offset.y  = -g_array_index(neighbours,Coordinates,i).y;

      for(j=0;j<input_bytes;j++) {
        offset.el = j;
        offset.offset = j
     - g_array_index(neighbours,Coordinates,i).x * mycorpus->depth 
     - g_array_index(neighbours,Coordinates,i).y * mycorpus->depth * mycorpus->width;
      g_array_append_val(Offsets,offset);
      }
    }

    offset.x  = 0;
    offset.y  = 0;
    offset.el = random_pos;
    offset.offset = random_pos;
    offset.weight = (int)(pow(2,16.0*parameters.random_weight));
    g_array_append_val(Offsets,offset);

    for(y=0;y<mydata->height;y++)
      for(x=0;x<mydata->width;x++)
        if (!bitmap_at(status,Status,x,y)[0].done)
     for(i=0;i<neighbours->len;i++) {
       int xx = x - g_array_index(neighbours,Coordinates,i).x;
       int yy = y - g_array_index(neighbours,Coordinates,i).y;

       bool use;
       if (!wrap_or_clip(&parameters,mydata,&xx,&yy))
         use = false; 
       else
         use = bitmap_at(status,Status,xx,yy)[0].done;
       
       if (use)
         bitmap_at(status,Status,x,y)[0].priority += (int)(weight_sum[i]*65536.0/sum); 
     }
  }

  g_free(weight_sum);


  qsort(& (g_array_index(Offsets,Offset,0)), Offsets->len, sizeof(Offset),
        offset_sorter);

  
  {
    Pixelel **items = g_new(Pixelel*,n_corpus_pixels);

    for(i=0,y=0;y<mycorpus->height;y++)
      for(x=0;x<mycorpus->width;x++)
        if (bitmap_at(mymask,Pixelel,x,y)[0])
          items[i++] = bitmap_at(mycorpus,Pixelel,x,y);

    tree = build_tree(items,n_corpus_pixels);  

    g_free(items);
  }

  /* ------------------------ */
  /* ------ Grow image ------ */
  /* ------------------------ */

  {
    Status *next;
    int n_pixels_done = 0;
    Pixelel* tmp_2 ;

    heap_become_heapish(heap);
  

   while(next = heap_pop_top(heap)) {
      Pixelel *result;
      int *weight;
      int pos = next-((Status*)status->data);
      int x = pos % status->width;
      int y = pos / status->width;

  
      
      g_array_set_size(Patterns,0);
      for(i=0;i<Offsets->len;i++) {
        Offset *o = &g_array_index(Offsets,Offset,i);
        int yy  = y + o->y;
        int xx  = x + o->x;
        int el = o->el;
      
        if (wrap_or_clip(&parameters,mydata,&xx,&yy))
     {
     Status *tmp = bitmap_at(status,Status,xx,yy) ;
     if ( tmp[0].use && (el >= input_bytes || tmp[0].done))
       {
            Pattern pattern;
            pattern.offset_no = i;
            pattern.value = bitmap_at(mydata,Pixelel,xx,yy)[el];
          g_array_append_val(Patterns,pattern);
            }
     }
      }
      
      bitmap_at(status,Status,x,y)[0].done = true;

      result = query_tree(tree);
      tmp_2 = bitmap_at(mydata,Pixelel,x,y) ;
      for(i=0;i<input_bytes;i++)
        tmp_2[i] = result[i];

      weight = (result - ((Pixelel*)mycorpus->data))/mycorpus->depth*weights->depth + weights->data;

      for(i=0;i<neighbours->len;i++) {
        int xx = x + g_array_index(neighbours,Coordinates,i).x;
        int yy = y + g_array_index(neighbours,Coordinates,i).y;

        if (wrap_or_clip(&parameters,mydata,&xx,&yy))
     {
     Status *tmp = bitmap_at(status,Status,xx,yy) ;
     if ( tmp[0].use && !tmp[0].done)
       {
            tmp[0].priority += weight[i];
            heap_promote(heap,tmp);
       }
          }
      }

      n_pixels_done++;
    } /* -----end while (!progress_monitor_aborted(monitor)...----- */
  } /* ------ end of Grow image ------ */

  node_free(tree);

  clean_up_tree_globals();

     for (i=0;i<mydata->width*mydata->height;i++){
      if (i%mydata->width == 0) printf("\n");
      printf("%i ",((Pixelel*)mydata->data)[i*4]);
   }

  bitmap_free(status);
  bitmap_free(weights);
  bitmap_free(mymask);
  bitmap_free(mycorpus);
  bitmap_free(mydata);

  heap_free(heap);

  g_array_free(neighbours,TRUE);
  
}

int main(int argC, char** argV){
   run();
}