std_forward_info.c

JPEG2000实现的源码

                   
{
    int out_neg, out_pos, in_taps, out_taps, synth_taps, i, k;
    float *sp, *op, *ip, *opp, val;
    
    out_neg = synth->neg_support*interp_factor + in->neg_support;
    out_pos = synth->pos_support*interp_factor + in->pos_support;
    out->neg_support = out_neg;
    out->pos_support = out_pos;
    out_taps = out_neg+out_pos+1;
    out->taps = (float *)
        local_realloc(out->taps,sizeof(float)*(size_t) out_taps);
    
    /* Perform interpolated convolution, as described above. */
    
    for (i=0; i < out_taps; i++)
        out->taps[i] = 0.0F; /* Ready to accumulate. */
    synth_taps = synth->neg_support + synth->pos_support + 1;
    in_taps = in->neg_support + in->pos_support + 1;
    op = NULL;
    for (sp=synth->taps,opp=out->taps,i=0;i<synth_taps;i++,opp+=interp_factor)
    {
        val = *(sp++);
        for (k=0, ip=in->taps, op=opp; k < in_taps; k++)
            *(op++) += val * *(ip++);
    }
    assert(op == (out->taps + out_taps));
}

/*****************************************************************************/
/* STATIC                     create_unit_waveform                           */
/*****************************************************************************/

static void
create_unit_waveform(std_tap_info_ptr taps)
{
    taps->neg_support = taps->pos_support = 0;
    taps->taps = (float *)
        local_malloc(INFO_MEM_KEY,sizeof(float));
    taps->taps[0] = 1.0F;
    taps->int_rdx = -1;
    taps->int_taps = NULL;
}

/*****************************************************************************/
/* STATIC                   compute_waveform_l2_norm                         */
/*****************************************************************************/

static float
compute_waveform_l2_norm(std_tap_info_ptr taps)
{
    double sum, val;
    int n, num_taps;
    
    num_taps = taps->neg_support + taps->pos_support + 1;
    for (sum=0.0, n=0; n < num_taps; n++)
    {
        val = taps->taps[n];
        sum += val*val;
    }
    return((float) sqrt(sum));
}

/*****************************************************************************/
/* STATIC                       get_nominal_gain                             */
/*****************************************************************************/

static float
get_nominal_gain(std_tap_info_ptr taps)

/* This function estimates the expansion in nominal range associated with
   a single one-dimensional stage in the Wavelet decomposition, by computing
   the absolute gain of the analysis kernel at DC and at Nyquist and
   returning the maximum of these two. */

{
    int n, num_taps;
    double dc, nyq, val;
    
    num_taps = taps->neg_support + taps->pos_support + 1;
    for (dc=nyq=0.0, n=0; n < num_taps; n++)
    {
        val = taps->taps[n];
        dc += val;
        nyq += (n&1)?val:(-val);
    }
    dc = fabs(dc);
    nyq = fabs(nyq);
    return((float)((dc>nyq)?dc:nyq));
}

/* SAIC General Decomp. Begin mods */
/*****************************************************************************/
/* STATIC                     compute_filter_gains                           */
/*****************************************************************************/

static void
compute_filter_gains(std_component_info_ptr comp_info)

/* This function traces the synthesis operations from each subband back to
   the original image, computing the expansion in nominal range (DC gain)
   and the expansion in squared error (wmse) and recording this information
   in the `nominal_gain' and `l2_norm' fields of each `band_info'
   structure. */

{
    std_level_info_ptr lev;
    std_band_info_ptr band;
    float nominal_gain;
    std_tap_info conv_h1,conv_h2,conv_h3,conv_v1,conv_v2,conv_v3,tmp_taps;
    int n,b,vert_hp_descent,horiz_hp_descent,max_hp_descent,d,orient;
    int level_depth,src_tree_depth,prev_vert_hp_descent,prev_horiz_hp_descent;
    int d_idx, d1, d2;

    /* Initialize the temporary synthesis waveform structures. */

    create_unit_waveform(&conv_h1);
    create_unit_waveform(&conv_h2);
    create_unit_waveform(&conv_h3);
    create_unit_waveform(&conv_v1);
    create_unit_waveform(&conv_v2);
    create_unit_waveform(&conv_v3);

    /* Now descend through the resolution levels, updating the base
       synthesis waveforms in `conv_h1' and `conv_v1' as we go to represent
       the LL band in each successive resolution level, and determining the
       waveforms for each of the level's high-pass bands as branches from this
       main path. */

    nominal_gain = 1.0F;
    for (level_depth=0, n=comp_info->num_levels; n >= 0; n--)
    {
        std_tap_info_ptr h_src, h_dst, v_src, v_dst;
        float hor_low_gain=0,vert_low_gain=0,hor_high_gain=0,vert_high_gain=0;

        lev = comp_info->levels + n;
        max_hp_descent = lev->max_hp_descent;
        if (n > 0)
        {
            hor_low_gain = get_nominal_gain(lev->hor_kernel.low);
            hor_high_gain = get_nominal_gain(lev->hor_kernel.high);
            vert_low_gain = get_nominal_gain(lev->vert_kernel.low);
            vert_high_gain = get_nominal_gain(lev->vert_kernel.high);
        }
        for (b=lev->min_band; b <= lev->max_band; b++)
        {
            band = lev->bands + b;
            if (!band->valid_band)
                continue;
            
            vert_hp_descent = band->vert_hp_descent;
            horiz_hp_descent = band->horiz_hp_descent;
            prev_vert_hp_descent = max_hp_descent;
            prev_horiz_hp_descent = max_hp_descent;
            
            h_dst = &conv_h1; v_dst = &conv_v1;
            band->nominal_gain = nominal_gain;
            for(src_tree_depth=level_depth,d_idx=0,
                d1=max_hp_descent-band->vert_hp_descent_chg[0],
                d2=max_hp_descent-band->horiz_hp_descent_chg[0]; 
                ((d1 >= max_hp_descent-vert_hp_descent) &&
                (d2 >= max_hp_descent-horiz_hp_descent)); 
                d1-=band->vert_hp_descent_chg[d_idx+1],
                d2-=band->horiz_hp_descent_chg[d_idx+1], 
                src_tree_depth++, d_idx++)
            {
                d = max_hp_descent - d_idx - 1;
                h_src = h_dst; v_src = v_dst;
                h_dst = (h_src==&conv_h2)?(&conv_h3):(&conv_h2);
                v_dst = (v_src==&conv_v2)?(&conv_v3):(&conv_v2);
                assert(b != LL_BAND);
                orient = (b>>(d+d)) & 3;
                switch (orient) {
                case LL_BAND:
                    if (d2 < prev_horiz_hp_descent) {
                        convolve_waveforms(lev->hor_kernel.low+1,h_src,h_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= hor_low_gain;
                    }
                    if (d1 < prev_vert_hp_descent) {
                        convolve_waveforms(lev->vert_kernel.low+1,v_src,v_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= vert_low_gain;
                    }
                    break;
                case HL_BAND:
                    if (d2 < prev_horiz_hp_descent) {
                        convolve_waveforms(lev->hor_kernel.high+1,h_src,h_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= hor_high_gain;
                    }
                    if (d1 < prev_vert_hp_descent) {
                        convolve_waveforms(lev->vert_kernel.low+1,v_src,v_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= vert_low_gain;
                    }
                    break;
                case LH_BAND:
                    if (d2 < prev_horiz_hp_descent) {
                        convolve_waveforms(lev->hor_kernel.low+1,h_src,h_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= hor_low_gain;
                    }
                    if (d1 < prev_vert_hp_descent) {
                        convolve_waveforms(lev->vert_kernel.high+1,v_src,v_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= vert_high_gain;
                    }
                    break;
                case HH_BAND:
                    if (d2 < prev_horiz_hp_descent) {
                        convolve_waveforms(lev->hor_kernel.high+1,h_src,h_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= hor_high_gain;
                    }
                    if (d1 < prev_vert_hp_descent) {
                        convolve_waveforms(lev->vert_kernel.high+1,v_src,v_dst,
                            1<<src_tree_depth);
                        band->nominal_gain *= vert_high_gain;
                    }
                    break;
                default: assert(0);
                }
                prev_vert_hp_descent = d1;
                prev_horiz_hp_descent = d2;
            }
            band->l2_norm =
                compute_waveform_l2_norm(h_dst) * compute_waveform_l2_norm(v_dst);
        }

        /* Now compute the next level's tree depth and base synthesis
           waveforms. */
        if (n == 0)
            break;
        h_src = &conv_h1; v_src = &conv_v1;
        h_dst = &conv_h2; v_dst = &conv_v2;
        convolve_waveforms(lev->hor_kernel.low+1,h_src,h_dst,1<<level_depth);
        convolve_waveforms(lev->vert_kernel.low+1,v_src,v_dst,1<<level_depth);
        tmp_taps = *h_src; *h_src = *h_dst; *h_dst = tmp_taps;
        tmp_taps = *v_src; *v_src = *v_dst; *v_dst = tmp_taps;
        level_depth++;
        nominal_gain *= hor_low_gain*vert_low_gain;
    }
    
    /* Cleanup. */
    
    destroy_taps(&conv_h1);
    destroy_taps(&conv_h2);
    destroy_taps(&conv_h3);
    destroy_taps(&conv_v1);
    destroy_taps(&conv_v2);
    destroy_taps(&conv_v3);
}
/* SAIC General Decomp. End mods */

/*****************************************************************************/
/* STATIC                   exponent_mantissa_to_step                        */
/*****************************************************************************/

static float
exponent_mantissa_to_step(int exponent, int mantissa)

/* Computes a relative or absolute step size as 2{-E} * (1+2^{-Mbits)*M),
   where E is the `exponent' value, M is the `mantissa' value and
   Mbits is the STEP_MANTISSA_BITS constant. */

{
    float step;
    
    assert((mantissa >= 0) && (mantissa < (1<<STEP_MANTISSA_BITS)));
    step = 1.0F + ((float) mantissa) / ((float)(1<<STEP_MANTISSA_BITS));
    while (exponent > 0)
    {
        step *= 0.5F;
        exponent--;
    }
    while (exponent < 0)
    {
        step *= 2.0F;
        exponent++;
    }
    return(step);
}

/*****************************************************************************/
/* STATIC                   step_to_exponent_mantissa                        */
/*****************************************************************************/

static void
step_to_exponent_mantissa(float step, int *exponent, int *mantissa)

/* Computes legal `exponent' and `mantissa' values which evaluate as closely
   as possible to `step' when supplied to `exponent_mantissa_to_step',
   as defined above. */

{
    int exp_val, mant_val;
    
    exp_val = 0;
    while (step >= 2.0F)
    {
        step *= 0.5F;
        exp_val--;
    }
    while (step < 1.0F)
    {
        step *= 2.0F;
        exp_val++;
    }
    mant_val = (int) floor((step-1.0)*((float)(1<<STEP_MANTISSA_BITS)) + 0.5);
    if (mant_val == (1<<STEP_MANTISSA_BITS))
    {
        mant_val = 0;