approx.c

linux下将各类格式图片转换工具

/*
 *  approx.c:		Approximation of range images with matching pursuit
 *
 *  Written by:		Ullrich Hafner
 *		
 *  This file is part of FIASCO (獸籸actal 獻籱age 獳籲d 玈籩quence 獵O籨ec)
 *  Copyright (C) 1994-2000 Ullrich Hafner <[e-mail]>
 */

/*
 *  $Date: 2000/06/14 20:50:51 $
 *  $Author: hafner $
 *  $Revision: 5.1 $
 *  $State: Exp $
 */

#include "config.h"

#include <math.h>

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

#include "cwfa.h"
#include "ip.h"
#include "rpf.h"
#include "domain-pool.h"
#include "misc.h"
#include "list.h"
#include "approx.h"
#include "coeff.h"
#include "wfalib.h"

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

			     local variables
  
*****************************************************************************/

typedef struct mp
{
   word_t exclude [MAXEDGES];
   word_t indices [MAXEDGES + 1];
   word_t into [MAXEDGES + 1];
   real_t weight [MAXEDGES];
   real_t matrix_bits;
   real_t weights_bits;
   real_t err;
   real_t costs;
} mp_t;

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

			     prototypes
  
*****************************************************************************/

static void
orthogonalize (unsigned index, unsigned n, unsigned level, real_t min_norm,
	       const word_t *domain_blocks, const coding_t *c);
static void 
matching_pursuit (mp_t *mp, bool_t full_search, real_t price,
		  unsigned max_edges, int y_state, const range_t *range,
		  const domain_pool_t *domain_pool, const coeff_t *coeff,
		  const wfa_t *wfa, const coding_t *c);

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

				public code
  
*****************************************************************************/

real_t 
approximate_range (real_t max_costs, real_t price, int max_edges,
		   int y_state, range_t *range, domain_pool_t *domain_pool,
		   coeff_t *coeff, const wfa_t *wfa, const coding_t *c)
/*
 *  Approximate image block 'range' by matching pursuit. This functions
 *  calls the matching pursuit algorithm several times (with different
 *  parameters) in order to find the best approximation. Refer to function
 *  'matching_pursuit()' for more details about parameters.
 *
 *  Return value:
 *	approximation costs
 */
{
   mp_t	  mp;
   bool_t success = NO;

   /*
    *  First approximation attempt: default matching pursuit algorithm.
    */
   mp.exclude [0] = NO_EDGE;
   matching_pursuit (&mp, c->options.full_search, price, max_edges,
		     y_state, range, domain_pool, coeff, wfa, c);

   /*
    *  Next approximation attempt: remove domain block mp->indices [0]
    *  from domain pool (vector with smallest costs) and run the
    *  matching pursuit again.
    */
   if (c->options.second_domain_block)
   {
      mp_t tmp_mp = mp;
      
      tmp_mp.exclude [0] = tmp_mp.indices [0];
      tmp_mp.exclude [1] = NO_EDGE;
	    
      matching_pursuit (&tmp_mp, c->options.full_search, price, max_edges,
			y_state, range, domain_pool, coeff, wfa, c);
      if (tmp_mp.costs < mp.costs)	/* success */
      {
	 success = YES;
	 mp      = tmp_mp;
      }
   }

   /*
    *  Next approximation attempt: check whether some coefficients have
    *  been quantized to zero. Vectors causing the underflow are
    *  removed from the domain pool and then the matching pursuit
    *  algorithm is run again (until underflow doesn't occur anymore).
    */
   if (c->options.check_for_underflow)
   {
      int  iteration = -1;
      mp_t tmp_mp    = mp;
      
      do
      {
	 int i;
 
	 iteration++;
	 tmp_mp.exclude [iteration] = NO_EDGE;
	 
	 for (i = 0; isdomain (tmp_mp.indices [i]); i++)
	    if (tmp_mp.weight [i] == 0)
	    {
	       tmp_mp.exclude [iteration] = tmp_mp.indices [i];
	       break;
	    }
      
	 if (isdomain (tmp_mp.exclude [iteration])) /* try again */
	 {
	    tmp_mp.exclude [iteration + 1] = NO_EDGE;
	    
	    matching_pursuit (&tmp_mp, c->options.full_search, price,
			      max_edges, y_state, range, domain_pool,
			      coeff, wfa, c);
	    if (tmp_mp.costs < mp.costs)	/* success */
	    {
	       success = YES;
	       mp      = tmp_mp;
	    }
	 }
      } while (isdomain (tmp_mp.exclude [iteration])
	       && iteration < MAXEDGES - 1);
   }

   /*
    *  Next approximation attempt: check whether some coefficients have
    *  been quantized to +/- max-value. Vectors causing the overflow are
    *  removed from the domain pool and then the matching pursuit
    *  algorithm is run again (until overflow doesn't occur anymore).
    */
   if (c->options.check_for_overflow)
   {
      int  iteration = -1;
      mp_t tmp_mp    = mp;
      
      do
      {
	 int i;
 
	 iteration++;
	 tmp_mp.exclude [iteration] = NO_EDGE;
	 
	 for (i = 0; isdomain (tmp_mp.indices [i]); i++)
	 {
	    rpf_t *rpf = tmp_mp.indices [i] ? coeff->rpf : coeff->dc_rpf;
	    
	    if (tmp_mp.weight [i] == btor (rtob (200, rpf), rpf)
		|| tmp_mp.weight [i] == btor (rtob (-200, rpf), rpf))
	    {
	       tmp_mp.exclude [iteration] = tmp_mp.indices [i];
	       break;
	    }
	 }
      
	 if (isdomain (tmp_mp.exclude [iteration])) /* try again */
	 {
	    tmp_mp.exclude [iteration + 1] = NO_EDGE;
	    
	    matching_pursuit (&tmp_mp, c->options.full_search, price,
			      max_edges, y_state, range, domain_pool,
			      coeff, wfa, c);
	    if (tmp_mp.costs < mp.costs)	/* success */
	    {
	       success = YES;
	       mp      = tmp_mp;
	    }
	 }
      } while (isdomain (tmp_mp.exclude [iteration])
	       && iteration < MAXEDGES - 1);
   }

   /*
    *  Finally, check whether the best approximation has costs
    *  smaller than 'max_costs'.
    */
   if (mp.costs < max_costs) 
   {
      int    edge;
      bool_t overflow  = NO;
      bool_t underflow = NO;
      int    new_index, old_index;

      new_index = 0;
      for (old_index = 0; isdomain (mp.indices [old_index]); old_index++)
	 if (mp.weight [old_index] != 0)
	 {
	    rpf_t *rpf = mp.indices [old_index] ? coeff->rpf : coeff->dc_rpf;
	    
	    if (mp.weight [old_index] == btor (rtob (200, rpf), rpf)
		|| mp.weight [old_index] == btor (rtob (-200, rpf), rpf))
	       overflow = YES;
	    
	    mp.indices [new_index] = mp.indices [old_index];
	    mp.into [new_index]    = mp.into [old_index];
	    mp.weight [new_index]  = mp.weight [old_index];
	    new_index++;
	 }
	 else
	    underflow = YES;
      
      mp.indices [new_index] = NO_EDGE;
      mp.into  [new_index]   = NO_EDGE;

      /*
       *  Update of probability models
       */
      {
	 word_t *domain_blocks = domain_pool->generate (range->level, y_state,
							wfa,
							domain_pool->model);
	 domain_pool->update (domain_blocks, mp.indices,
			      range->level, y_state, wfa, domain_pool->model);
	 coeff->update (mp.weight, mp.into, range->level, coeff);
	 
	 Free (domain_blocks);
      }
      
      for (edge = 0; isedge (mp.indices [edge]); edge++)
      {
	 range->into [edge]   = mp.into [edge];
	 range->weight [edge] = mp.weight [edge];
      }
      range->into [edge]  = NO_EDGE;
      range->matrix_bits  = mp.matrix_bits;
      range->weights_bits = mp.weights_bits;
      range->err          = mp.err;
   }
   else
   {
      range->into [0] = NO_EDGE;
      mp.costs	      = MAXCOSTS;
   }
   
   return mp.costs;
}

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

			     local variables
  
*****************************************************************************/

static real_t norm_ortho_vector [MAXSTATES];     
/*
 *  Square-norm of the i-th vector of the orthogonal basis (OB)
 *  ||o_i||^2; i = 0, ... ,n
 */
static real_t ip_image_ortho_vector [MAXEDGES];
/* 
 *  Inner product between the i-th vector of the OB and the given range: 
 *  <b, o_i>; i = 0, ... ,n 
 */
static real_t ip_domain_ortho_vector [MAXSTATES][MAXEDGES];
/* 
 *  Inner product between the i-th vector of the OB and the image of domain j: 
 *  <s_j, o_i>; j = 0, ... , wfa->states; i = 0, ... ,n, 
 */
static real_t rem_denominator [MAXSTATES];     
static real_t rem_numerator [MAXSTATES];
/*
 *  At step n of the orthogonalization the comparitive value
 *  (numerator_i / denominator_i):= <b, o_n>^2 / ||o_n|| ,
 *  is computed for every domain i,
 *  where o_n := s_i - \sum(k = 0, ... , n-1) {(<s_i, o_k> / ||o_k||^2) o_k}
 *  To avoid computing the same values over and over again,
 *  the constant (remaining) parts of every domain are
 *  stored in 'rem_numerator' and 'rem_denominator' separately
 */
static bool_t used [MAXSTATES];    
/*
 *  Shows whether a domain image was already used in a
 *  linear combination (YES) or not (NO)
 */

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

				private code
  
*****************************************************************************/

static void 
matching_pursuit (mp_t *mp, bool_t full_search, real_t price,
		  unsigned max_edges, int y_state, const range_t *range,
		  const domain_pool_t *domain_pool, const coeff_t *coeff,
		  const wfa_t *wfa, const coding_t *c)
/*
 *  Find an approximation of the current 'range' with a linear
 *  combination of vectors of the 'domain_pool'. The linear
 *  combination is generated step by step with the matching pursuit
 *  algorithm.  If flag 'full_search' is set then compute complete set
 *  of linear combinations with n = {0, ..., 'max_edges'} vectors and
 *  return the best one. Otherwise abort the computation as soon as
 *  costs (LC (n + 1)) exceed costs ( LC (n)) and return the
 *  sub-optimal solution.  'price' is the langrange multiplier
 *  weighting rate and distortion.  'band' is the current color band
 *  and 'y_state' the corresponding state in the Y component at same
 *  pixel position.  'domain_pool' gives the set of available vectors,
 *  and 'coeff' the model for the linear factors. The number of
 *  elements in the linear combination is limited by 'max_edges'. In
 *  'mp', vectors may be specified which should be excluded during the
 *  approximation.
 *  
 *  No return value.
 *
 *  Side effects:
 *	vectors, factors, rate, distortion and costs are stored in 'mp'
 */
{
   unsigned	 n;			/* current vector of the OB */
   int		 index;			/* best fitting domain image */
   unsigned	 domain;		/* counter */
   real_t	 norm;			/* norm of range image */
   real_t	 additional_bits;	/* bits for mc, nd, and tree */
   word_t	*domain_blocks;		/* current set of domain images */
   const real_t  min_norm = 2e-3;	/* lower bound of norm */
   unsigned 	 best_n   = 0;