component_mix.c

JPEG2000实现的源码

    int i, j, k=0;
    
    if (tile_kernel->third_d_kernel_user_defined) {
      std_user_defined_ptr user;
      
      user = std_user_defined_summary(&(tile_kernel->third_d_kernel));
      /* Don't send entire user structure, just what is needed for this
     particular filter. */
      data_size = 1 + 1 + 2 * user->num_steps;
      for (i=0; i<user->num_steps; i++) {
    data_size += 2 * user->step_supports[i];
      }
      data = (char *)local_malloc(MIXING_MEM_KEY,data_size*sizeof(char));
      data[k++] = tile_kernel->third_d_kernel_user_defined;
      data[k++] = user->num_steps;
      for (i=0; i<user->num_steps; i++) {
    data[k++] = user->step_supports[i];
    data[k++] = user->downshifts[i];
      }
      for (i=0; i<user->num_steps; i++) {
    for (j=0; j<user->step_supports[i]; j++) {
      data[k++] = (user->taps[i][j] & 0x0000FF00) >> 8;
      data[k++] = (user->taps[i][j] & 0x000000FF);
    }
      }
      local_free(user);
    }
    else {
      data_size = 1 + strlen(tile_kernel->third_d_identifier) + 1;
      data = (char *) local_malloc(MIXING_MEM_KEY, data_size * sizeof(char));
      data[0] = tile_kernel->third_d_kernel_user_defined;
      memcpy(data+1,tile_kernel->third_d_identifier, 
         strlen(tile_kernel->third_d_identifier) + 1);
    }
    self->stream->set_main_cmetag(self->stream, 
                  (std_ushort) self->next_xform_marker, 
                  data_size,
                  (std_byte *) data);
    local_free(data);
  }
  
  return(kernel_type);
}

/*****************************************************************************/
/* STATIC                 get_third_d_convolution_taps                       */
/*****************************************************************************/

static void
get_third_d_convolution_taps(the_component_mix_ref self,
                 float **low_taps, float **high_taps)
{
  component_tile_kernel_ptr tile_kernel;
  std_tap_info_ptr low = 0, high = 0;
  int kernel_type = 0;
  
  if (self->current_tile_idx >= self->num_tiles)
    local_error("Attempting to access 3-D DWT information beyond "
        "the last image tile!");
  tile_kernel = self->tile_kernels + self->current_tile_idx;
  
  low = tile_kernel->third_d_kernel.low;
  high = tile_kernel->third_d_kernel.high;
  kernel_type = tile_kernel->third_d_kernel.kernel_type;
  
  if (kernel_type != INFO__CONVOLUTION)
    local_error("Attempting to retrieve convolution taps for a "
        "Wavelet kernel which is to be implemented by lifting "
        "rather than convolution!!");
  if (low_taps != NULL)
    *low_taps = low[0].taps + low[0].neg_support;
  if (high_taps != NULL)
    *high_taps = high[0].taps + high[0].neg_support;
}

/*****************************************************************************/
/* STATIC                   get_third_d_int_steps                            */
/*****************************************************************************/

static int
get_third_d_int_steps(the_component_mix_ref self, int neg_support[],
              int pos_support[], int *taps[], int downshift[])
{
  component_tile_kernel_ptr tile_kernel;
  std_lifting_info_ptr info = 0;
  int n, kernel_type = 0;
  
  if (self->current_tile_idx >= self->num_tiles)
    local_error("Attempting to access 3-D DWT information beyond "
        "the last image tile!");
  tile_kernel = self->tile_kernels + self->current_tile_idx;
  
  info = &(tile_kernel->third_d_kernel.lifting);
  kernel_type = tile_kernel->third_d_kernel.kernel_type;
  
  if (kernel_type != INFO__INT_LIFTING)
    local_error("Attempting to retrieve reversible lifting steps "
        "for a Wavelet kernel which is to be implemented using "
        "either non-reversible lifting steps or convolution!");
  assert(info->num_steps <= INFO__MAX_LIFTING_STEPS);
  for (n=0; n < info->num_steps; n++)
    {
      neg_support[n] = info->steps[n].neg_support;
      pos_support[n] = info->steps[n].pos_support;
      taps[n] = info->steps[n].int_taps + info->steps[n].neg_support;
      downshift[n] = info->steps[n].int_rdx;
    }
  return(info->num_steps);
}

/*****************************************************************************/
/* STATIC                         apply_int_step                             */
/*****************************************************************************/

static void
apply_int_step(int the_stage,int neg_support,int pos_support,int *taps,
           int downshift,ifc_int *src,ifc_int *dst,int src_cols,
           int dst_cols,int src_offset,int src_on_left_boundary,
           int src_on_right_boundary,int skip_first,int skip_last)
{
  int even_taps, n;
  int left_ext, right_ext;
  std_int sum, offset;
  ifc_int *sp,*dp;
  
  assert(src_cols && dst_cols);
  even_taps = (pos_support+neg_support+1)&1;
  if (even_taps){ /* Extension for even tap filters. */
    if (src_offset)
      local_error("Cannot start image, tile or frame with an odd indexed "
          "sample when using even length filters!  Use odd length "
          "filters!");
    left_ext = symmetric_extension_type(1,(src_cols+dst_cols)&1,
                    the_stage&1,0);
    right_ext = symmetric_extension_type(1,(src_cols+dst_cols)&1,
                     the_stage&1,1);
  }
  else{ /* Extension for odd tap filters. */
    pos_support-=src_offset;
    neg_support+=src_offset;
    taps += src_offset;
    left_ext = (src_on_left_boundary)?0:1;
    right_ext = (src_on_right_boundary)?0:1;
  }
  sp = src; dp = sp + src_cols - 1;
  /* extension to the left*/
  for(n=1;n<=neg_support;n++){
    /* Half point extension*/
    if(left_ext){
      assert( n  <= src_cols );  
      /* TJF: above was formerly assert( n  <= src_cols ) */
      sp[-n] = sp[n-1]*left_ext;
    }
    /* whole point extesnion*/
    else{
      assert( n <= src_cols);
      sp[-n] = sp[n];
    }
  }
  /* Extension to the right*/
  for(n=1;n<=pos_support+dst_cols-src_cols;n++){
    /* 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 extesnion*/
    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                       perform_1d_int_lifting                       */
/*****************************************************************************/

static void
perform_1d_int_lifting(lifting_info_ref self, int buffer_length)
     /* Performs 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=0;n<self->num_steps;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;
    apply_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                     check_kernel_symmetry                          */
/*****************************************************************************/

static void
check_kernel_symmetry(int neg_support,int pos_support,float *taps,
              int even_taps,int high_pass)
     /* Checks that the supplied kernel has the relevant symmetry (even or odd,
    as identified by the `even_taps' argument).  Note that we can be
    quite sure that the kernel has been normalized to have a unit
    DC or Nyquist gain, as appropriate, which allows us to establish
    reliable tollerances.  If the test fails, an appropriate error
    message is generated. */
{
  double tolerance;
  int failed;
  int n;
  
  tolerance = 0.0001;
  failed = 0;
  if(even_taps){
    if ((neg_support-high_pass) != (pos_support+high_pass-1))
      failed = 1;
    if(high_pass && !failed){
      for (n=1; n <= neg_support; n++)
    if (fabs(taps[-n]+taps[n-1]) > tolerance)
      failed = 1;
    }
    if ((!high_pass) && !failed){
      for (n=1; n <= pos_support; n++)
    if (fabs(taps[n]-taps[-n+1]) > tolerance)
      failed = 1;
    }
  }
  else{
    if (neg_support != pos_support)
      failed = 1;
    if (!failed)
      for (n=1; n <= pos_support; n++)
    if (fabs(taps[n]-taps[-n]) > tolerance)
      failed = 1;
  }
  if (failed)
    local_error("Wavelet kernels may be implemented by convolution "
        "only if they satisfy the symmetry requirements imposed "
        "by the boundary extension policy.  All other kernels "
        "must be implemented by lifting -- i.e. the identifier "
        "supplied via the `-Fkernels' argument must be preceded "
        "with a capital `L'.");
}

/*****************************************************************************/
/*                              perform_1d_analysis                          */
/*****************************************************************************/

static void
perform_1d_analysis(conv_info_ref base, int length)
     /* This function symmetrically extends the supplied input buffer as necessary
    and writes the two sets of subband samples into the two buffers pointed to
    by the entries of the `branch_bufs' array. */
{
  float *sp, *dp;
  int n, b;
  
  if (length == 1)
    { /* Generally speaking, the transforms only work for signals of length
     2 or more. */
      *(base->branch_bufs[0]) = base->in_buf[0];
      return;
    }
  
  /* Perform symmetric extension.  This is almost completely trivial, but
     we must beware of the possibility that the entire row is smaller than
     the kernel support, even down to a length of 1.  As a consequence,
     the best way to go about applying the extension is to extend each
     boundary one sample at a time. */
  
  sp = base->in_buf; dp = sp + length - 1;
  if(base->even_taps){
    for(n=1;n<=base->extend;n++){
      sp[-n] = sp[n-1];
      dp[n] = dp[-n+1];
    }
  }
  else{ /* The usual case of odd-length filters is the most general and also
       the simplest. */
    for (n=1;n<=base->extend;n++){
      sp[-n] = sp[n];
      dp[n] = dp[-n];
    }
  }
  
  /* Apply the two kernels one at a time. */
  
  for(b=0;b<2;b++){
    int neg, pos, k;
    float *taps, sum;
    
    dp = base->branch_bufs[b];