path: root/converter/other/fiasco/codec/subdivide.c

/*
 *  subdivide.c:        Recursive subdivision of range images
 *
 *  Written by:         Ullrich Hafner
 *
 *  This file is part of FIASCO (Fractal Image And Sequence COdec)
 *  Copyright (C) 1994-2000 Ullrich Hafner
 */

/*
 *  $Date: 2000/07/15 17:59:31 $
 *  $Author: hafner $
 *  $Revision: 5.4 $
 *  $State: Exp $
 */

#include "config.h"

#include <string.h>

#include "pm_c_util.h"

#include "types.h"
#include "macros.h"
#include "error.h"

#include "image.h"
#include "cwfa.h"
#include "approx.h"
#include "ip.h"
#include "bintree.h"
#include "control.h"
#include "prediction.h"
#include "domain-pool.h"
#include "mwfa.h"
#include "misc.h"
#include "subdivide.h"
#include "list.h"
#include "coeff.h"
#include "wfalib.h"

/*****************************************************************************

                                prototypes

*****************************************************************************/

static void
init_new_state (bool_t auxiliary_state, bool_t delta, range_t *range,
                const range_t *child, const int *y_state,
                wfa_t *wfa, coding_t *c);
static void
init_range (range_t *range, const image_t *image, unsigned band,
            const wfa_t *wfa, coding_t *c);

/*****************************************************************************

                                public code

*****************************************************************************/

real_t
subdivide (real_t max_costs, unsigned band, int y_state, range_t *range,
           wfa_t *wfa, coding_t *c, bool_t prediction, bool_t delta)
/*
 *  Subdivide the current 'range' recursively and decide whether
 *  a linear combination, a recursive subdivision, or a prediction is
 *  the best choice of approximation.
 *  'band' is the current color band, 'y_state' is the corresponding
 *  state of the Y color component (color image compression only).
 *  If 'prediction' is TRUE then also test motion compensation or
 *  nondeterministic approximation.
 *  If 'delta' is TRUE then current range is already predicted.
 *
 *  Return value:
 *      costs of the best approximation or MAXCOSTS if costs exceed 'max_costs'
 *
 *  Side effects:
 *      'range' factors and costs of linear combination are modified
 *      'wfa'   new transitions and prediction coefficients are added
 *      'c'     pixels and inner products are updated
 */
{
   real_t    subdivide_costs;        /* Costs arising from approx. the current
                                       range with two children */
   real_t    lincomb_costs;          /* Costs arising from approx. the current
                                       range with a linear combination */
   int       new_y_state [MAXLABELS];   /* Corresponding state of Y */
   real_t    price;                     /* Approximation costs multiplier */
   bool_t    try_mc;                    /* YES: try MC prediction */
   bool_t    try_nd;                    /* YES: try ND prediction */
   unsigned  states;                    /* Number of states before the
                                           recursive subdivision starts */
   void     *domain_model;              /* copy of domain pool model */
   void     *d_domain_model;            /* copy of delta domain pool model */
   void     *lc_domain_model;           /* copy of domain pool model */
   void     *lc_d_domain_model;         /* copy of delta domain pool model */
   void     *coeff_model;               /* copy of coefficients model */
   void     *d_coeff_model;             /* copy of delta coefficients model */
   void     *lc_coeff_model;            /* copy of coefficients model */
   void     *lc_d_coeff_model;          /* copy of delta coefficients model */
   tree_t    tree_model;                /* copy of tree model */
   tree_t    p_tree_model;              /* copy of pred. tree model */
   range_t   lrange;                    /* range of lin. comb. approx. */
   range_t   rrange;                    /* range of recursive approx. */
   range_t   child [MAXLABELS];         /* new children of the current range */
   static unsigned percent = 0;         /* status of progress meter */

   if (wfa->wfainfo->level == range->level)
      percent = 0;

   range->into [0] = NO_EDGE;           /* default approximation: empty */
   range->tree     = RANGE;

   if (range->level < 3)                /* Don't process small ranges */
      return MAXCOSTS;

   /*
    *  If image permutation (tiling) is performed and the tiling level
    *  is reached then get coordinates of the new block.
    */
   if (c->tiling->exponent
       && range->level == wfa->wfainfo->level - c->tiling->exponent)
   {
      unsigned width, height;           /* size of range (dummies)*/

      if (c->tiling->vorder [range->global_address] < 0)
         return 0;                      /* nothing to do */
      else
         locate_subimage (wfa->wfainfo->level, range->level,
                          c->tiling->vorder [range->global_address],
                          &range->x, &range->y, &width, &height);
   }

   if (range->x >= c->mt->original->width ||
       range->y >= c->mt->original->height)
      return 0;                         /* range is not visible */

   /*
    *  Check whether prediction is allowed or not
    *  mc == motion compensation, nd == nondeterminism
    */
   try_mc = (prediction && c->mt->frame_type != I_FRAME
             && range->level >= wfa->wfainfo->p_min_level
             && range->level <= wfa->wfainfo->p_max_level
             && (range->x + width_of_level (range->level)
                 <= c->mt->original->width)
             && (range->y + height_of_level (range->level)
                 <= c->mt->original->height));

   try_nd = (prediction && c->mt->frame_type == I_FRAME
             && range->level >= wfa->wfainfo->p_min_level
             && range->level <= wfa->wfainfo->p_max_level);

   if (try_mc)
      clear_norms_table (range->level, wfa->wfainfo, c->mt);


   /*
    *  Check if current range must be initialized. I.e. range pixels must
    *  be copied from entire image to bintree pixel buffer. Moreover,
    *  all inner products tables must be initialized.
    */
   if (range->level == c->options.lc_max_level)
      init_range (range, c->mt->original, band, wfa, c);

   price = c->price;
   if (band != Y)
      price *= c->options.chroma_decrease; /* less quality for chroma bands */

   /*
    *  Compute children of corresponding state in Y band
    */
   if (band != Y)                       /* Cb and Cr bands only */
   {
      unsigned label;

      for (label = 0; label < MAXLABELS; label++)
         if (ischild (y_state))
            new_y_state [label] = wfa->tree [y_state][label];
         else
            new_y_state [label] = RANGE;
   }
   else
      new_y_state [0] = new_y_state [1] = RANGE;

   /*
    *  Store contents of all models that may get modified during recursion
    */
   domain_model   = c->domain_pool->model_duplicate (c->domain_pool->model);
   d_domain_model = c->d_domain_pool->model_duplicate (c->d_domain_pool->model);
   coeff_model    = c->coeff->model_duplicate (c->coeff, c->coeff->model);
   d_coeff_model  = c->d_coeff->model_duplicate (c->d_coeff, c->d_coeff->model);
   tree_model     = c->tree;
   p_tree_model   = c->p_tree;
   states         = wfa->states;

   /*
    *  First alternative of range approximation:
    *  Compute costs of linear combination.
    */
   if (range->level <= c->options.lc_max_level) /* range is small enough */
   {
      lrange                 = *range;
      lrange.tree            = RANGE;
      lrange.tree_bits       = tree_bits (LEAF, lrange.level, &c->tree);
      lrange.matrix_bits     = 0;
      lrange.weights_bits    = 0;
      lrange.mv_tree_bits    = try_mc ? 1 : 0; /* mc allowed but not used */
      lrange.mv_coord_bits   = 0;
      lrange.nd_tree_bits    = 0;
      lrange.nd_weights_bits = 0;
      lrange.prediction      = NO;

      lincomb_costs
         = approximate_range (max_costs, price, c->options.max_elements,
                              y_state, &lrange,
                              (delta ? c->d_domain_pool : c->domain_pool),
                              (delta ? c->d_coeff : c->coeff), wfa, c);
   }
   else
      lincomb_costs = MAXCOSTS;

   /*
    *  Store contents of models that have been modified
    *  by approximate_range () above ...
    */
   lc_domain_model   = c->domain_pool->model;
   lc_d_domain_model = c->d_domain_pool->model;
   lc_coeff_model    = c->coeff->model;
   lc_d_coeff_model  = c->d_coeff->model;
   /*
    *  ... and restore them with values before lc
    */
   c->domain_pool->model   = c->domain_pool->model_duplicate (domain_model);
   c->d_domain_pool->model = c->d_domain_pool->model_duplicate (d_domain_model);
   c->coeff->model         = c->coeff->model_duplicate (c->coeff, coeff_model);
   c->d_coeff->model       = c->d_coeff->model_duplicate (c->d_coeff,
                                                          d_coeff_model);

   /*
    *  Second alternative of range approximation:
    *  Compute costs of recursive subdivision.
    */
   if (range->level > c->options.lc_min_level) /* range is large enough */
   {
      unsigned label;

      memset (&child [0], 0, 2 * sizeof (range_t)); /* initialize children */

      /*
       *  Initialize a new range for recursive approximation
       */
      rrange                 = *range;
      rrange.tree_bits       = tree_bits (CHILD, rrange.level, &c->tree);
      rrange.matrix_bits     = 0;
      rrange.weights_bits    = 0;
      rrange.err             = 0;
      rrange.mv_tree_bits    = try_mc ? 1 : 0;  /* mc allowed but not used */
      rrange.mv_coord_bits   = 0;
      rrange.nd_tree_bits    = try_nd ?
                               tree_bits (CHILD, lrange.level, &c->p_tree): 0;
      rrange.nd_weights_bits = 0;
      rrange.prediction      = NO;

      /*
       *  Initialize the cost function and subdivide the current range.
       *  Every child is approximated by a recursive call of subdivide()
       */
      subdivide_costs = (rrange.tree_bits + rrange.weights_bits
                         + rrange.matrix_bits + rrange.mv_tree_bits
                         + rrange.mv_coord_bits + rrange.nd_tree_bits
                         + rrange.nd_weights_bits) * price;

      for (label = 0; label < MAXLABELS; label++)
      {
         real_t remaining_costs;        /* upper limit for next recursion */

         child[label].image          = rrange.image * MAXLABELS + label + 1;
         child[label].address        = rrange.address * MAXLABELS + label;
         child[label].global_address = rrange.global_address * MAXLABELS
                                       + label;
         child[label].level          = rrange.level - 1;
         child[label].x = rrange.level & 1
                          ? rrange.x
                          : (rrange.x
                             + label * width_of_level (rrange.level - 1));
         child[label].y = rrange.level & 1
                          ? (rrange.y
                             + label * height_of_level (rrange.level - 1))
                          : rrange.y;

         /*
          *  If necessary compute the inner products of the new states
          *  (generated during the recursive approximation of child [0])
          */
         if (label && rrange.level <= c->options.lc_max_level)
            compute_ip_images_state (child[label].image, child[label].address,
                                     child[label].level, 1, states, wfa, c);
         /*
          *  Call subdivide() for both children.
          *  Abort the recursion if 'subdivide_costs' exceed 'lincomb_costs'
          *  or 'max_costs'.
          */
         remaining_costs = MIN(lincomb_costs, max_costs) - subdivide_costs;

         if (remaining_costs > 0)       /* still a way for improvement */
         {
            subdivide_costs += subdivide (remaining_costs, band,
                                          new_y_state [label], &child [label],
                                          wfa, c, prediction, delta);
         }
         else if (try_mc && child[label].level >= wfa->wfainfo->p_min_level)
         {
            fill_norms_table (child[label].x, child[label].y,
                              child[label].level, wfa->wfainfo, c->mt);
         }

         if (try_mc)
            update_norms_table (rrange.level, wfa->wfainfo, c->mt);

         /*
          *  Update of progress meter
          */
         if (c->options.progress_meter != FIASCO_PROGRESS_NONE)
         {
            if (c->options.progress_meter == FIASCO_PROGRESS_PERCENT)
            {
               unsigned new_percent;    /* new status of progress meter */

               new_percent = (child[label].global_address + 1) * 100.0
                             / (1 << (wfa->wfainfo->level - child[label].level));
               if (new_percent > percent)
               {
                  percent = new_percent;
                  info ("%3d%%  \r", percent);
               }
            }
            else if (c->options.progress_meter == FIASCO_PROGRESS_BAR)
            {
               unsigned new_percent;    /* new status of progress meter */

               new_percent = (child[label].global_address + 1) * 50.0
                             / (1 << (wfa->wfainfo->level
                                      - child[label].level));
               for (; new_percent > percent; percent++)
               {
                  info ("#");
               }
            }
         }

         /*
          *  If costs of subdivision exceed costs of linear combination
          *  then abort recursion.
          */
         if (subdivide_costs >= MIN(lincomb_costs, max_costs))
         {
            subdivide_costs = MAXCOSTS;
            break;
         }
         rrange.err             += child [label].err;
         rrange.tree_bits       += child [label].tree_bits;
         rrange.matrix_bits     += child [label].matrix_bits;
         rrange.weights_bits    += child [label].weights_bits;
         rrange.mv_tree_bits    += child [label].mv_tree_bits;
         rrange.mv_coord_bits   += child [label].mv_coord_bits;
         rrange.nd_weights_bits += child [label].nd_weights_bits;
         rrange.nd_tree_bits    += child [label].nd_tree_bits;

         tree_update (ischild (child [label].tree) ? CHILD : LEAF,
                      child [label].level, &c->tree);
         tree_update (child [label].prediction ? LEAF : CHILD,
                      child [label].level, &c->p_tree);
      }
   }
   else
      subdivide_costs = MAXCOSTS;

   /*
    *  Third alternative of range approximation:
    *  Predict range via motion compensation or nondeterminism and
    *  approximate delta image.
    */
   if (try_mc || try_nd)                /* try prediction */
   {
       real_t prediction_costs; /* Costs arising from approx. the current
                                   range with prediction */

       prediction_costs
           = predict_range (MIN(MIN(lincomb_costs, subdivide_costs),
                                max_costs),
                            price, range, wfa, c, band, y_state, states,
                            &tree_model, &p_tree_model, domain_model,
                            d_domain_model, coeff_model, d_coeff_model);
       if (prediction_costs < MAXCOSTS) /* prediction has smallest costs */
       {
           c->domain_pool->model_free (domain_model);
           c->d_domain_pool->model_free (d_domain_model);
           c->domain_pool->model_free (lc_domain_model);
           c->d_domain_pool->model_free (lc_d_domain_model);
           c->coeff->model_free (coeff_model);
           c->d_coeff->model_free (d_coeff_model);
           c->coeff->model_free (lc_coeff_model);
           c->d_coeff->model_free (lc_d_coeff_model);

           return prediction_costs;
       }
   }

   if (lincomb_costs >= MAXCOSTS && subdivide_costs >= MAXCOSTS)
   {
      /*
       *  Return MAXCOSTS if neither a linear combination nor a recursive
       *  subdivision yield costs less than 'max_costs'
       */
      c->domain_pool->model_free (c->domain_pool->model);
      c->d_domain_pool->model_free (c->d_domain_pool->model);
      c->domain_pool->model_free (lc_domain_model);
      c->d_domain_pool->model_free (lc_d_domain_model);

      c->coeff->model_free (c->coeff->model);
      c->d_coeff->model_free (c->d_coeff->model);
      c->coeff->model_free (lc_coeff_model);
      c->d_coeff->model_free (lc_d_coeff_model);

      c->domain_pool->model   = domain_model;
      c->d_domain_pool->model = d_domain_model;
      c->coeff->model         = coeff_model;
      c->d_coeff->model       = d_coeff_model;
      c->tree                 = tree_model;
      c->p_tree               = p_tree_model;

      if (wfa->states != states)
         remove_states (states, wfa);

      return MAXCOSTS;
   }
   else if (lincomb_costs < subdivide_costs)
   {
      /*
       *  Use the linear combination: The factors of the linear combination
       *  are stored already in 'range', so revert the probability models
       *  only.
       */
      c->domain_pool->model_free (c->domain_pool->model);
      c->d_domain_pool->model_free (c->d_domain_pool->model);
      c->domain_pool->model_free (domain_model);
      c->d_domain_pool->model_free (d_domain_model);

      c->coeff->model_free (c->coeff->model);
      c->d_coeff->model_free (c->d_coeff->model);
      c->coeff->model_free (coeff_model);
      c->d_coeff->model_free (d_coeff_model);

      c->domain_pool->model   = lc_domain_model;
      c->d_domain_pool->model = lc_d_domain_model;
      c->coeff->model         = lc_coeff_model;
      c->d_coeff->model       = lc_d_coeff_model;
      c->tree                 = tree_model;
      c->p_tree               = p_tree_model;

      *range = lrange;

      if (wfa->states != states)
         remove_states (states, wfa);

      return lincomb_costs;
   }
   else
   {
      /*
       *  Use the recursive subdivision: Generate a new state with transitions
       *  given in child[].
       *  Don't use state in linear combinations in any of the following cases:
       *  - if color component is Cb or Cr
       *  - if level of state > tiling level
       *  - if state is (partially) outside image geometry
       */
      if (band > Y
          || (c->tiling->exponent
              && rrange.level > wfa->wfainfo->level - c->tiling->exponent)
          || (range->x + width_of_level (range->level)
              > c->mt->original->width)
          || (range->y + height_of_level (range->level)
              > c->mt->original->height))
         init_new_state (YES, delta, &rrange, child, new_y_state, wfa, c);
      else
         init_new_state (NO, delta, &rrange, child, new_y_state, wfa, c);

      *range = rrange;

      c->domain_pool->model_free (domain_model);
      c->d_domain_pool->model_free (d_domain_model);
      c->domain_pool->model_free (lc_domain_model);
      c->d_domain_pool->model_free (lc_d_domain_model);
      c->coeff->model_free (coeff_model);
      c->d_coeff->model_free (d_coeff_model);
      c->coeff->model_free (lc_coeff_model);
      c->d_coeff->model_free (lc_d_coeff_model);

      return subdivide_costs;
   }
}

void
cut_to_bintree (real_t *dst, const word_t *src,
                unsigned src_width, unsigned src_height,
                unsigned x0, unsigned y0, unsigned width, unsigned height)
/*
 *  Cut region ('x0', 'y0', 'width', 'height') of the pixel array 'src'.
 *  Size of image is given by 'src_width' x 'src_height'.
 *  'dst' pixels are converted to bintree order and real format.
 *
 *  No return value.
 *
 *  Side effects:
 *      'dst []' is filled with corresponding region.
 */
{
   const unsigned mask01      = 0x555555; /* binary ...010101010101 */
   const unsigned mask10      = 0xaaaaaa; /* binary ...101010101010 */
   const unsigned mask01plus1 = mask01 + 1; /* binary ...010101010110 */
   const unsigned mask10plus1 = mask10 + 1; /* binary ...101010101011 */
   unsigned       x, y;                 /* pixel coordinates */
   unsigned       xmask, ymask;         /* address conversion */

   if (width != height && width != (height >> 1))
      error ("Bintree cutting requires special type of images.");

   ymask = 0;
   for (y = y0; y < y0 + height; y++, ymask = (ymask + mask10plus1) & mask01)
   {
      xmask = 0;
      for (x = x0; x < x0 + width; x++, xmask = (xmask + mask01plus1) & mask10)
      {
         if (y >= src_height || x >= src_width)
            dst [xmask | ymask] = 0;
         else
            dst [xmask | ymask] = src [y * src_width + x] / 16;
      }
   }
}

/*****************************************************************************

                                private code

*****************************************************************************/

static void
init_new_state (bool_t auxiliary_state, bool_t delta, range_t *range,
                const range_t *child, const int *y_state,
                wfa_t *wfa, coding_t *c)
/*
 *  Initializes a new state with all parameters needed for the encoding step.
 *  If flag 'auxiliary_state' is set then don't insert state into domain pools.
 *  If flag 'delta' is set then state represents a delta image (prediction via
 *  nondeterminism or motion compensation).
 *  'range' the current range image,
 *   'child []' the left and right children of 'range'.
 *
 *  No return value.
 *
 *  Side effects:
 *      New state is appended to 'wfa' (and also its inner products and images
 *      are computed and stored in 'c')
 */
{
   unsigned label;
   bool_t   state_is_domain = NO;

   if (!auxiliary_state)
   {
      if (!delta || c->options.delta_domains)
         state_is_domain = c->domain_pool->append (wfa->states, range->level,
                                                   wfa, c->domain_pool->model);
      if (delta || c->options.normal_domains)
         state_is_domain = c->d_domain_pool->append (wfa->states, range->level,
                                                     wfa,
                                                     c->d_domain_pool->model)
                           || state_is_domain;
   }
   else
      state_is_domain = NO;

   range->into [0] = NO_EDGE;
   range->tree     = wfa->states;

   for (label = 0; label < MAXLABELS; label++)
   {
      wfa->tree [wfa->states][label]       = child [label].tree;
      wfa->y_state [wfa->states][label]    = y_state [label];
      wfa->mv_tree [wfa->states][label]    = child [label].mv;
      wfa->x [wfa->states][label]          = child [label].x;
      wfa->y [wfa->states][label]          = child [label].y;
      wfa->prediction [wfa->states][label] = child [label].prediction;

      append_transitions (wfa->states, label, child [label].weight,
                          child [label].into, wfa);
   }
   wfa->delta_state [wfa->states] = delta;

   if (range->err < 0)
      warning ("Negative image norm: %f, %f", child [0].err, child [1].err);

/*    state_is_domain = YES; */

   append_state (!state_is_domain,
                 compute_final_distribution (wfa->states, wfa),
                 range->level, wfa, c);
}

static void
init_range (range_t *range, const image_t *image, unsigned band,
            const wfa_t *wfa, coding_t *c)
/*
 *  Read a new 'range' of the image 'image_name' (current color component
 *  is 'band') and compute the new inner product arrays.
 *
 *  No return value.
 *
 *  Side effects:
 *      'c->pixels' are filled with pixel values of image block
 *      'c->ip_images_state' are computed with respect to new image block
 *      'range->address' and 'range->image' are initialized with zero
 */
{
   unsigned state;

   /*
    *  Clear already computed products
    */
   for (state = 0; state < wfa->states; state++)
      if (need_image (state, wfa))
         memset (c->ip_images_state[state], 0,
                 size_of_tree (c->products_level) * sizeof(real_t));

   cut_to_bintree (c->pixels, image->pixels [band],
                   image->width, image->height,
                   range->x, range->y, width_of_level (range->level),
                   height_of_level (range->level));

   range->address = range->image = 0;
   compute_ip_images_state (0, 0, range->level, 1, 0, wfa, c);
}