hplx_hor.c

JPEG2000 EBCOT算法源码

/*****************************************************************************/
/* Copyright 1998, Hewlett-Packard Company                                   */
/* All rights reserved                                                       */
/* Author: Christos Chrysafis                                                */
/* Version: V2.1                                                             */
/* Last Revised: 10/2/98                                                     */
/*****************************************************************************/

#include <local_heap.h>
#include <memory.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <line_block_ifc.h>
#include <hplx_hor.h>
#include <hplx_coeffs.h>
#include <hplx_def.h>




/*******************************************************************************/
/*  inner product two pixels at a time */
/******************************************************************************/
void two_at_once_h(int plus, int *flag,int i,float m,float *a,float *b,float *c)
{
	if(!m) return;

	if(plus)
	{
		if(*flag)
			while(i--)    *a++  = m * ( *(b+=2) + *(c+=2));
			else
				while(i--)    *a++ += m * ( *(b+=2) + *(c+=2));
	}
	else
	{
		if(*flag)
			while(i--)    *a++  = m * ( *(b+=2) - *(c+=2));
			else
				while(i--)    *a++ += m * ( *(b+=2) - *(c+=2));
	}
	*flag=0;
}
/*******************************************************************************/
/*  inner product one pixel at a time */
/*******************************************************************************/
void one_at_once_h(int *flag,int i,float m,float *a,float *b)
{
	if(!m) return;
	
		if(*flag)
			while(i--)    *a++  = m * *(b+=2) ;
			else
				while(i--)    *a++ += m *  *(b+=2) ;
	
	*flag=0;
}

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

/************************************************************************/
static  int 
inner_prod_fixed( int * a,  int  b,
                 short int  *c, int  k)
{       
  while( b-- )
    k+=((int) *(a++)) * ((int)*(c++));
  return(k);
}

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

static float
inner_prod(float* a, float *b, float *c,float k)
{       
  while( a != b )
    k+=*(a++)**(c++);
  return(k);
}

/************************************************************************/
/* Set the size for the one dimensional line 
* to be decomposed.
*
*
*  One object can be used for different decompositions
*
*/
/************************************************************************/

static void
decimator__set_size(decimator_ref self, int t)
{
  self->Length=min(t,self->xdim);
}      

/*****************************************************************************/
/* STATIC                  compute_shift_factor                              */
/*****************************************************************************/

static int
compute_shift_factor(float *l_taps, int l_num, float *h_taps, int h_num,
                     int precision, float max_scale)
{
  float val, max_val;
  int i;
  
  max_val = 0.0F;
  for (i=0; i < l_num; i++)
  {
    val = l_taps[i];
    if (val > max_val)
      max_val = val;
    else if ((-val) > max_val)
      max_val = -val;
  }
  for (i=0; i < h_num; i++)
  {
    val = h_taps[i];
    if (val > max_val)
      max_val = val;
    else if ((-val) > max_val)
      max_val = -val;
  }
  max_val *= max_scale;
  max_val *= 2.0F; /* all scaled coefficients will fit within this dynamic
                   range.  Want 2^R >= `max_val'*2^S where R is the
  precision and S is the shift we want to find. */
  val = ((float)(1<<precision)) / max_val;
  i=0;
  while (val >= 1.0F)
  {
    val *= 0.5F;
    i++;
  }
  i--;
  assert(((1<<i)*max_val <= (float)(1<<precision)) &&
    ((1<<i)*max_val > (float)(1<<(precision-1))));
  return(i);
}

/************************************************************************/
/*
*  Destructor for decimator object
*
*/
/************************************************************************/

static int
decimator__terminate(decimator_ref self )
{
  if(self->fixed)
  {
    local_free(self->buffer_I-=self->V_max);
    local_free(self->BufFast);
    local_free(self->VL_Low_I);
    local_free(self->VH_Low_I);
  }
  else
  {
    local_free(self->buffer-=self->V_max);
    local_free(self->VL_Low);
    local_free(self->VH_Low);
  }
  
  local_free(self);
  return(0);  
}

/************************************************************************/
/* Initialization routine for class decimator
*
*/
/************************************************************************/

static void
  decimator__initialize(decimator_ref self,      
                        int  X,
                        filter_info_ref info,
                        int c_levs, int component_idx)
{       
  int i,PLUSMINUS1;
  float *l_taps, *h_taps;
  int h_neg, h_pos, l_neg, l_pos, Lmin, Lbeyond, Hmin, Hbeyond;
  char *built_in_id;
  
  self->DefaultFB = 0;
  self->component_idx = component_idx;
  
  l_taps = info->get_taps(info,0,component_idx,0x0F&c_levs,&l_neg,&l_pos,NULL);
  h_taps = info->get_taps(info,FILTER__HIGH,component_idx,0x0F&c_levs,
                          &h_neg,&h_pos,&built_in_id);
  if ((built_in_id != NULL) && (strcmp(built_in_id,"0") == 0))
    self->DefaultFB = 1;
  
  self->EvenFilter = (l_neg+l_pos) & 1;
  Lmin = -l_pos; Lbeyond = l_neg+1;
  Hmin = -h_pos; Hbeyond = h_neg+1;
  if (self->EvenFilter)
  { Hmin--; Hbeyond--; }
  
  self->VL_min = Lmin;
  self->VL_max = Lbeyond;
  self->VH_min = Hmin;
  self->VH_max = Hbeyond;
  self->V_min = min( self->VL_min, self->VH_min );
  self->V_max = max( self->VL_max, self->VH_max );
  PLUSMINUS1     = self->EvenFilter?(-1):1;          
  

  self->VL_Low = (float *)
    local_malloc(sizeof(float)*(size_t)(Lbeyond-Lmin));
  self->VH_Low = (float *)
    local_malloc(sizeof(float)*(size_t)(Hbeyond-Hmin));
  self->VL_High = self->VL_Low + (Lbeyond - Lmin);
  self->VH_High = self->VH_Low + (Hbeyond - Hmin);
  for (i=0; i < (Lbeyond - Lmin); i++)
    self->VL_High[-i-1] = l_taps[i];
  for (i=0; i < (Hbeyond - Hmin); i++)
    self->VH_High[-i-1] = h_taps[i];
  
  for(i=0;i<self->VH_max-self->VH_min;i++)
    self->VH_Low[i] *=PLUSMINUS1;
  
  self->xdim=self->Length=X;
         
  self->buffer =
    ((float*)
    local_malloc((self->Length-self->V_min+self->V_max)*sizeof(float)))
    + self->V_max;
}

/************************************************************************/
/*
*  decomposes an input line, 
*
*/
/************************************************************************/

static void 
decimator__decimate(decimator_ref self, float * V)
{

  int j,i,evenF,effectiveL,flag;
  float *SindxL,*SindxH,*Sindex1,*OindxH,qu=0;
  float *downL,*downH;
  int  shift_factorL, shift_factorH,_Strange,OddSignal;  

  downL=self->VL_Low;
  downH=self->VH_Low;
  evenF=self->EvenFilter;
 
  effectiveL=self->Length;
  if( (effectiveL&1) && evenF)
  {effectiveL--; _Strange=1;} else _Strange=0;
  
  /***********************************************
  ** This is the main part of the signal
  ***********************************************/
  
  OddSignal = effectiveL & 1;
  shift_factorL =  - self->VL_max - 1;
  shift_factorH =  - self->VH_max - evenF ;

  if (_Strange)
	  qu =(float)( (*(self->buffer+effectiveL))); 
  
  /***********************************************
  ** Extension to the left of the signal
  ***********************************************/
  
  for(j=-1,i=-(int)(evenF);j>-self->V_max;j--)
	  self->buffer[j]=self->buffer[(++i)%effectiveL];
  
  /***********************************************
  ** Extension to the right of the signal
  ***********************************************/
  
  for ( j = effectiveL ,
	  i = ( ( effectiveL << 1 ) - !evenF ) ;
  j < effectiveL - self->V_min; j++ )  
	  self->buffer[j]=self->buffer[(--i)%effectiveL];
  
  /********************************
  Filtering for the Low pass
  Filtering for the High pass
  ********************************/
  
      SindxL        = self->buffer + shift_factorL ;
      OindxH       = V + ((effectiveL+1+_Strange)>>1);
      SindxH        = self->buffer + shift_factorH ;
      
 if(!self->DefaultFB)
   {   
	  /* Low Pass*/
	  j=(effectiveL>>1)+OddSignal;
	  flag=1;
      i=self->VL_max-self->VL_min;
	  Sindex1=SindxL+i-1;
	  i>>=1;
	  while(i--)
          two_at_once_h(1, &flag,j, *(downL++) ,V, SindxL++,Sindex1--);

      if(!evenF)   /*Oddlength Filter */
          one_at_once_h(&flag,j,*(downL++),V,SindxL++);
	  
	  V+=j;
      if(_Strange) *(V++)=qu;  
      /*  High Pass */  

	  j=(effectiveL>>1);
	  flag=1;
      i=self->VH_max-self->VH_min;
	  Sindex1=SindxH+i-1;
	  i>>=1;

	  while(i--)
        two_at_once_h(!evenF, &flag,j, *(downH++) ,V, SindxH++,Sindex1--);

       if(!evenF)   /*Oddlength Filter */
          one_at_once_h(&flag,j,*(downH++),V,SindxH++);
      
      } 
      else
      {  
        for(i=0;i<(effectiveL>>1);i++)
        {
          SindxL+=2;
          *(V++) =     
            (float)(ha4) * (SindxL[0]+SindxL[8])+
            (float)(ha3) * (SindxL[1]+SindxL[7])+
            (float)(ha2) * (SindxL[2]+SindxL[6])+
            (float)(ha1) * (SindxL[3]+SindxL[5])+
            (float)(ha0) *  SindxL[4] ;
          *(OindxH++)=   
            (float)(ga3) * (SindxL[2]+SindxL[8])+
            (float)(ga2) * (SindxL[3]+SindxL[7])+
            (float)(ga1) * (SindxL[4]+SindxL[6])+ 
            (float)(ga0) *  SindxL[5];
        }		 
  
        if(OddSignal)
        {
          SindxL+=2;
          *(V++) =     
            (float)(ha4) * (SindxL[0]+SindxL[8])+
            (float)(ha3) * (SindxL[1]+SindxL[7])+
            (float)(ha2) * (SindxL[2]+SindxL[6])+
            (float)(ha1) * (SindxL[3]+SindxL[5])+
            (float)(ha0) *  SindxL[4] ;
        }
		 if(_Strange) *(V)=qu;   
      }
}



/************************************************************************/
static void
interpolator__set_size(interpolator_ref self, int t)
{
  self->Length=min(t,self->xdim);
  
}


/************************************************************************/
/*
*  initialization routine for class interpolator
*
*   c_levs           levels of decomposition performed
*/
/************************************************************************/

static void
  interpolator__initialize(interpolator_ref self, 
                           int  X, filter_info_ref info,int c_levs,
                           int component_idx)
{
  int i,PLUSMINUS1;
  float *Lfilt,*Hfilt;
  float *l_taps, *h_taps; 
  int h_neg, h_pos, l_neg, l_pos, Lmin, Lbeyond, Hmin, Hbeyond;
  char *built_in_id;
  
  self->DefaultFB = 0;
  self->component_idx = component_idx;
 
  l_taps = info->get_taps(info,FILTER__SYNTH,component_idx,0x0F & c_levs,
                          &l_neg,&l_pos,NULL);
  h_taps = info->get_taps(info,FILTER__HIGH | FILTER__SYNTH,component_idx,
                          0x0F & c_levs,&h_neg,&h_pos,&built_in_id);
  if (strcmp(built_in_id,"0") == 0)
    self->DefaultFB = 1;
  
  self->EvenFilter = (l_neg+l_pos) & 1;
  Lmin = -l_neg; Lbeyond = l_pos+1;
  Hmin = -h_neg; Hbeyond = h_pos+1;
  if (self->EvenFilter)
  { Lmin--; Lbeyond--; }
  
  self->VH_min = Hmin;
  self->VH_max = Hbeyond;
  self->VL_min = Lmin;
  self->VL_max = Lbeyond;
  self->V_min = min( self->VL_min, self->VH_min );
  self->V_max = max( self->VL_max, self->VH_max );
  
  self->_data = (float*)
    local_malloc(sizeof(float)*(X));
  
  self->xdim = self->Length = X;
 
  PLUSMINUS1 =  self->EvenFilter?(-1):1; 
  
  i= max(max(self->VL_max,abs(self->VL_min)),
    max(self->VH_max,abs(self->VH_min)));
  
 
  self->buffer =
    ((float*)
    local_malloc((self->Length-self->V_min + self->V_max)*sizeof(float)))
    + self->V_max;

  Lfilt =
    ((float*)
    local_malloc((self->V_max-self->V_min+1)*sizeof(float)))
    - self->V_min;
  Hfilt =
    ((float*)
    local_malloc((self->V_max-self->V_min+1)*sizeof(float)))
    - self->V_min;
  
  for(i = self->V_min; i <self->V_max ; i++)
    Lfilt[i]=Hfilt[i]=0;
  
  l_taps -= Lmin;
  h_taps -= Hmin; /* Get centres. */
  
  for(i=0;i<self->V_max;i+=2)
  {
    if(i-self->EvenFilter<self->VH_max)
      Hfilt[-i]=h_taps[i-self->EvenFilter];
    else
      Hfilt[-i]=(float)(0.0);
    
    if(i<self->VL_max)
      Lfilt[-i-self->EvenFilter]=l_taps[i];
    else