/*
 *  ip.c:		Computation of inner products
 *
 *  Written by:		Ullrich Hafner
 *
 *  This file is part of FIASCO (Fractal Image And Sequence COdec)
 *  Copyright (C) 1994-2000 Ullrich Hafner
 */

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

#include "config.h"

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

#include "cwfa.h"
#include "control.h"
#include "ip.h"

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

				prototypes

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

static real_t
standard_ip_image_state (unsigned address, unsigned level, unsigned domain,
			 const coding_t *c);
static real_t
standard_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
			 const coding_t *c);

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

				public code

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

real_t
get_ip_image_state (unsigned image, unsigned address, unsigned level,
		    unsigned domain, const coding_t *c)
/*
 *  Return value:
 *	Inner product between 'image' ('address') and
 *      'domain' at given 'level'
 */
{
   if (level <= c->options.images_level)
   {
      /*
       *  Compute the inner product in the standard way by multiplying
       *  the pixel-values of the given domain and range image.
       */
      return standard_ip_image_state (address, level, domain, c);
   }
   else
   {
      /*
       *  Use the already computed inner products stored in 'ip_images_states'
       */
      return c->ip_images_state [domain][image];
   }
}

void
compute_ip_images_state (unsigned image, unsigned address, unsigned level,
			 unsigned n, unsigned from,
			 const wfa_t *wfa, coding_t *c)
/*
 *  Compute the inner products between all states
 *  'from', ... , 'wfa->max_states' and the range images 'image'
 *  (and children) up to given level.
 *
 *  No return value.
 *
 *  Side effects:
 *	inner product tables 'c->ip_images_states' are updated
 */
{
   if (level > c->options.images_level)
   {
      unsigned state, label;

      if (level > c->options.images_level + 1)	/* recursive computation */
	 compute_ip_images_state (MAXLABELS * image + 1, address * MAXLABELS,
				  level - 1, MAXLABELS * n, from, wfa, c);

      /*
       *  Compute inner product <f, Phi_i>
       */
      for (label = 0; label < MAXLABELS; label++)
	 for (state = from; state < wfa->states; state++)
	    if (need_image (state, wfa))
	    {
	       unsigned  edge, count;
	       int     	 domain;
	       real_t 	*dst, *src;

	       if (ischild (domain = wfa->tree [state][label]))
	       {
		  if (level > c->options.images_level + 1)
		  {
		     dst = c->ip_images_state [state] + image;
		     src = c->ip_images_state [domain]
			   + image * MAXLABELS + label + 1;
		     for (count = n; count; count--, src += MAXLABELS)
			*dst++ += *src;
		  }
		  else
		  {
		     unsigned newadr = address * MAXLABELS + label;

		     dst = c->ip_images_state [state] + image;

		     for (count = n; count; count--, newadr += MAXLABELS)
			*dst++ += standard_ip_image_state (newadr, level - 1,
							   domain, c);
		  }
	       }
	       for (edge = 0; isedge (domain = wfa->into [state][label][edge]);
		    edge++)
	       {
		  real_t weight = wfa->weight [state][label][edge];

		  if (level > c->options.images_level + 1)
		  {
		     dst = c->ip_images_state [state] + image;
		     src = c->ip_images_state [domain]
			   + image * MAXLABELS + label + 1;
		     for (count = n; count; count--, src += MAXLABELS)
			*dst++ += *src * weight;
		  }
		  else
		  {
		     unsigned newadr = address * MAXLABELS + label;

		     dst = c->ip_images_state [state] + image;

		     for (count = n; count; count--, newadr += MAXLABELS)
			*dst++ += weight *
				  standard_ip_image_state (newadr, level - 1,
							   domain, c);
		  }
	       }
	    }
   }
}

real_t
get_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
		    const coding_t *c)
/*
 *  Return value:
 *	Inner product between 'domain1' and 'domain2' at given 'level'.
 */
{
   if (level <= c->options.images_level)
   {
      /*
       *  Compute the inner product in the standard way by multiplying
       *  the pixel-values of both state-images
       */
      return standard_ip_state_state (domain1, domain2, level, c);
   }
   else
   {
      /*
       *  Use already computed inner products stored in 'ip_images_states'
       */
      if (domain2 < domain1)
	 return c->ip_states_state [domain1][level][domain2];
      else
	 return c->ip_states_state [domain2][level][domain1];
   }
}

void
compute_ip_states_state (unsigned from, unsigned to,
			 const wfa_t *wfa, coding_t *c)
/*
 *  Computes the inner products between the current state 'state1' and the
 *  old states 0,...,'state1'-1
 *
 *  No return value.
 *
 *  Side effects:
 *	inner product tables 'c->ip_states_state' are computed.
 */
{
   unsigned level;
   unsigned state1, state2;

   /*
    *  Compute inner product <Phi_state1, Phi_state2>
    */

   for (level = c->options.images_level + 1;
	level <= c->options.lc_max_level; level++)
      for (state1 = from; state1 <= to; state1++)
	 for (state2 = 0; state2 <= state1; state2++)
	    if (need_image (state2, wfa))
	    {
	       unsigned	label;
	       real_t	ip = 0;

	       for (label = 0; label < MAXLABELS; label++)
	       {
		  int	   domain1, domain2;
		  unsigned edge1, edge2;
		  real_t   sum, weight2;

		  if (ischild (domain1 = wfa->tree [state1][label]))
		  {
		     sum = 0;
		     if (ischild (domain2 = wfa->tree [state2][label]))
			sum = get_ip_state_state (domain1, domain2,
						  level - 1, c);

		     for (edge2 = 0;
			  isedge (domain2 = wfa->into [state2][label][edge2]);
			  edge2++)
		     {
			weight2 = wfa->weight [state2][label][edge2];
			sum += weight2 * get_ip_state_state (domain1, domain2,
							     level - 1, c);
		     }
		     ip += sum;
		  }
		  for (edge1 = 0;
		       isedge (domain1 = wfa->into [state1][label][edge1]);
		       edge1++)
		  {
		     real_t weight1 = wfa->weight [state1][label][edge1];

		     sum = 0;
		     if (ischild (domain2 = wfa->tree [state2][label]))
			sum = get_ip_state_state (domain1, domain2,
						  level - 1, c);

		     for (edge2 = 0;
			  isedge (domain2 = wfa->into [state2][label][edge2]);
			  edge2++)
		     {
			weight2 = wfa->weight [state2][label][edge2];
			sum += weight2 * get_ip_state_state (domain1, domain2,
							     level - 1, c);
		     }
		     ip += weight1 * sum;
		  }
	       }
	       c->ip_states_state [state1][level][state2] = ip;
	    }
}

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

				private code

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

static real_t
standard_ip_image_state (unsigned address, unsigned level, unsigned domain,
			 const coding_t *c)
/*
 *  Returns the inner product between the subimage 'address' and the
 *  state image 'domain' at given 'level'.  The stored state images
 *  and the image tree are used to compute the inner product in the
 *  standard way by multiplying the corresponding pixel values.
 *
 *  Return value:
 *	computed inner product
 */
{
   unsigned i;
   real_t   ip = 0, *imageptr, *stateptr;

   if (level > c->options.images_level)
      error ("We cannot interpret a Level %d image.", level);

   imageptr = &c->pixels [address * size_of_level (level)];

   stateptr = c->images_of_state [domain] + address_of_level (level);

   for (i = size_of_level (level); i; i--)
      ip += *imageptr++ * *stateptr++;

   return ip;
}

static real_t
standard_ip_state_state (unsigned domain1, unsigned domain2, unsigned level,
			 const coding_t *c)
/*
 *  Returns the inner product between the subimage 'address' and the
 *  state image 'state' at given 'level'.  The stored state images are
 *  used to compute the inner product in the standard way by
 *  multiplying the corresponding pixel values.
 *
 *  Return value:
 *	computed inner product
 */
{
   unsigned i;
   real_t   ip = 0, *state1ptr, *state2ptr;

   if (level > c->options.images_level)
      error ("We cannot interpret and image with Level %d.", level);

   state1ptr = c->images_of_state [domain1] + address_of_level (level);
   state2ptr = c->images_of_state [domain2] + address_of_level (level);

   for (i = size_of_level (level); i; i--)
      ip += *state1ptr++ * *state2ptr++;

   return ip;
}