component_demix.c

JPEG2000实现的源码

    {
      /* Half point extension*/
      if(right_ext)
        {
          assert( src_cols >=n-1);
          /* TJF: above was formerly assert( src_cols >=n) */
          dp[n] = dp[-n+1]*right_ext;
        }
      /* whole point extension*/
      else
        {
          assert(  src_cols - 1 >= n);
          dp[n] = dp[-n];
        }
    }

  /* Apply lifting step.  Note that `taps' points to "centre" of support. */

  offset = (1<<downshift) >> 1;
  if (skip_first)
    { dst_cols--; dst++; src++; }
  if (skip_last)
    { dst_cols--; }
  for (; dst_cols > 0; dst_cols--, src++, dst++)
    {
      for (sum=0, n=-neg_support; n <= pos_support; n++)
        sum += src[n] * taps[n];
      *dst -= (ifc_int)((sum+offset) >> downshift);
    }
}


/*****************************************************************************/
/* STATIC                        undo_1d_int_lifting                         */
/*****************************************************************************/

static void
  undo_1d_int_lifting(lifting_info_ref self, int buffer_length)

 /* Undoes the sequence of integer lifting steps on the
    two branch buffers in the `self->branch_buf_int' array. */

{
  ifc_int *work_bufs[2];
  int b, n, first_col_idx;
  int special_first, special_last;
  int left_boundary_odd, right_boundary_odd;
  int subseq_cols[2]; /* even and odd sub-sequence lengths. */

  /* Set up the `work_bufs' array first. */

  for (b=0; b < 2; b++)
    {
      work_bufs[b] = self->branch_buf_int[b];
    }

  /* Now for the main processing loop. */

  first_col_idx = self->first_idx;

  special_first = special_last = 0;
  if (buffer_length == 1)
    return;

  left_boundary_odd = first_col_idx & 1;
  right_boundary_odd = (first_col_idx + buffer_length - 1) & 1;
  subseq_cols[0] = (buffer_length+1-left_boundary_odd) >> 1;
  subseq_cols[1] = (buffer_length+left_boundary_odd) >> 1;
  for (n=self->num_steps-1; n >= 0; n--)
    {
      int src_seq, dst_seq;
      int src_offset; /* First src index - first dst index */

      src_seq = n & 1; dst_seq = (n+1) & 1;
      src_offset = (src_seq&1)?(-left_boundary_odd):left_boundary_odd;
      undo_int_step(n,self->neg_supports[n],self->pos_supports[n],
                    self->int_taps[n],self->int_downshifts[n],
                    work_bufs[src_seq],work_bufs[dst_seq],
                    subseq_cols[src_seq],subseq_cols[dst_seq],
                    src_offset,left_boundary_odd^dst_seq,
                    right_boundary_odd^dst_seq,
                    special_first&src_seq,special_last&src_seq);
    }
}

/*****************************************************************************/
/* STATIC                    create_polyphase_kernel                         */
/*****************************************************************************/

static void
  create_polyphase_kernel(hplx_polyphase_kernel_ptr kernel,
                          int input_odd, int output_odd,
                          int neg_support, int pos_support, float *taps)

 /* This function fills out the contents of the polyphase kernel structure
    referenced by `kernel', by reversing and sub-sampling the impulse
    response supplied via `taps', which may be indexed as `taps'[k] where
    k ranges from -`neg_support' to +`pos_support'.  The `input_odd'
    argument identifies whether this impulse response is notionally
    centred on even (low-pass) or odd (high-pass) samples, while the
    `output_odd' argument identifies whether the polyphase kernel is to
    generate the even or odd sub-sequence of the reconstructed signal. */

{
  int displacement, n;

  if (output_odd)
    output_odd = 1; /* Make sure. */
  if (input_odd)
    input_odd = 1; /* Make sure. */
  displacement = output_odd - input_odd;
  kernel->pos_support = (displacement + neg_support)>>1;
  kernel->neg_support = (pos_support - displacement)>>1;
  if ((kernel->neg_support < 0) || (kernel->pos_support < 0))
    return;
  kernel->taps = (float *)
    local_malloc(MIXING_MEM_KEY,
                 sizeof(float)*(kernel->neg_support+kernel->pos_support+1));
  kernel->taps += kernel->neg_support;
  for (n=-kernel->neg_support; n <= kernel->pos_support; n++)
    kernel->taps[n] = taps[displacement - 2*n];
}


/*****************************************************************************/
/* STATIC                      perform_1d_synthesis                          */
/*****************************************************************************/

static void
  perform_1d_synthesis(conv_info_ref base, int length)

 /* This function symmetrically extends the two input subband
    buffers in the `branch_bufs' array as necessary and writes the
    synthesized sequence into the buffer pointed to by `out_buf'. */

{
  int b, n;
  float *sp, *dp;

  /* First, apply the symmetric extension policy described in "ifc.h". */

  if (base->even_taps)
    {
      for (b=0; b < 2; b++)
        {
          int odd_right, band_cols;
          float sign_flip;

          sign_flip = (b==1)?-1.0F:1.0F;
          band_cols = (length+1-b)>>1;
          sp = base->branch_bufs[b];
          dp = sp + band_cols - 1;
          odd_right = (length & 1);
          if (odd_right && b)
            *(dp++) = 0.0F;
          for (n=1; n <= base->extend; n++)
            { sp[-n] = sign_flip * sp[n-1];
              dp[n] =  sign_flip * dp[-n+1-odd_right]; }
        }
    }
  else
    { /* Odd length filters with symmetric extension.  This is the nice
         easy case, which also is most flexible. */
      for (b=0; b < 2; b++)
        {
          int odd_left, odd_right, band_cols;

          odd_left = 1 - (b);
          odd_right = (length + b)  & 1;
          band_cols = (length+1-b)>>1;
          sp = base->branch_bufs[b];
          dp = sp + band_cols - 1;
          if (odd_left) /* Interleave right extensions before left. */
            for (n=1; n <= base->extend; n++)
              { dp[n] = dp[-n+1-odd_right];
                sp[-n] = sp[n-1+odd_left]; }
          else /* Interleave left extensions before right. */
            for (n=1; n <= base->extend; n++)
              { sp[-n] = sp[n-1+odd_left];
                dp[n] = dp[-n+1-odd_right]; }
        }
    }

  /* Now apply each of the four polyphase synthesis kernels. */

  for (b=0; b < 2; b++)
    { /* Generate the b'th output sub-sequence: b=0 for even; b=1 for odd. */
      hplx_polyphase_kernel_ptr kernel;
      float *taps, sum;
      int subseq_cols, neg, pos, k;

      if (b == 0)
        kernel = base->branch_to_even;
      else
        kernel = base->branch_to_odd;
      subseq_cols = (length+1-b)>>1;

      /* Apply kernel from low-pass band to get initial output samples. */

      sp = base->branch_bufs[0];
      neg = kernel->neg_support;
      pos = kernel->pos_support;
      taps = kernel->taps;
      for (dp=base->out_buf+(b),
           n=subseq_cols; n > 0; n--, dp+=2, sp++)
        {
          for (sum=0.0F, k=-neg; k <= pos; k++)
            sum += taps[k] * sp[k];
          *dp = sum;
        }

      /* Apply kernel from high-pass band and add to initial output samples. */

      kernel++;
      sp = base->branch_bufs[1];
      neg = kernel->neg_support;
      pos = kernel->pos_support;
      taps = kernel->taps;
      for (dp=base->out_buf+(b),
           n=subseq_cols; n > 0; n--, dp+=2, sp++)
        {
          for (sum=0.0F, k=-neg; k <= pos; k++)
            sum += taps[k] * sp[k];
          *dp += sum;
        }
    }
}
/* UofA End */

/* Kodak/SAIC Linear Transform Begin */
/*****************************************************************************/
/* STATIC                        invert_lin_xform                            */
/*****************************************************************************/

static void invert_lin_xform(ifc_int **components,
                             float **T_rev, int num_components, int width)
                             
{
    float *components_in, sum;
    int i, j;
    
    components_in = (float *) local_malloc(MIXING_MEM_KEY, 
        num_components*sizeof(float));
    for(width--;width>=0;width--){ /* Bug fix by Fred Wheeler */
        for(i=0;i<num_components;i++){
            components_in[i] = (float) components[i][width];
        }
        for(i=0;i<num_components;i++){
            sum=0;
            for(j=0;j<num_components;j++){
                sum+=components_in[j]*T_rev[i][j]; /* Bug fix by Fred Wheeler */
            }
            components[i][width] = (ifc_int) (sum + 0.5);
        }
    }
    local_free(components_in);
}
/* Kodak/SAIC Linear Transform End */

/* Begin Aerospace MCT mods (TSW) */
/*****************************************************************************/
/* STATIC                       invert_seg_lin_xform                         */
/*****************************************************************************/

static void
invert_seg_lin_xform(ifc_int **buffers,int num_components,int width,
             component_xfms_ref transforms,
             collection_ref seg_collection)
{
  /* This module allows for the number of input and output components
     to be different. It also handles specification of null matrices
     in the collection segments.*/
  
  float *interm_components,**interm_ptr,**xform_mat,*offset_vec;
  int i,j,k,l,n_in,n_out;
  ifc_int **comp_ptr;
  segment seg;
  
  interm_components = (float *) local_malloc(MIXING_MEM_KEY, 
                         num_components*width*sizeof(float));
  comp_ptr = (ifc_int **) local_malloc(MIXING_MEM_KEY, 
                       num_components*sizeof(ifc_int *));
  interm_ptr = (float **) local_malloc(MIXING_MEM_KEY, 
                       num_components*sizeof(float *));
  
  /* do all of the decorrelation transforms, putting the result 
     into the intermediate component buffers */
  
  for(i=0;i<seg_collection->n_comp_segments;i++){
    seg = seg_collection->comp_segment[i];
    n_in = seg.n_input_components;
    
    for(j=0;j<n_in;j++)
      comp_ptr[j]=buffers[seg.input_components[j]];
    n_out=seg.n_output_components;
    for(j=0;j<n_out;j++)
      interm_ptr[j]=interm_components+seg.output_components[j]*width;
    xform_mat=transforms->decorrel_xfms[seg.xfm_matrix].matrix_values;
    offset_vec=transforms->decorrel_offs[seg.off_vector].vector_values;
    
    /* It is probably most efficient to deal with null transform matrices
       here by using 4 different pieces of code for the 4 possible cases,
       otherwise much time will be spent allocating and deallocating 
       memory */
    
    if((seg.xfm_matrix>0)&&(seg.off_vector>0)){
      for(j=0;j<width;j++){
    for(k=0;k<n_out;k++){
      *interm_ptr[k] = 0;
      for(l=0;l<n_in;l++)
        *interm_ptr[k] +=
          (((float) *comp_ptr[l])*xform_mat[k][l]);
      *(interm_ptr[k]++) += offset_vec[k];
    }
    for(l=0;l<n_in;l++)
      comp_ptr[l]++;
      }
    }
    else if((seg.xfm_matrix>0)&&(seg.off_vector==0)){
      for(j=0;j<width;j++){
    for(k=0;k<n_out;k++){
      *interm_ptr[k] = 0;
      for(l=0;l<n_in;l++)
        *interm_ptr[k] +=