001/*
002 * $RCSfile: AnWTFilterFloatLift9x7.java,v $
003 * $Revision: 1.1 $
004 * $Date: 2005/02/11 05:02:29 $
005 * $State: Exp $
006 *
007 * Class:                   AnWTFilterFloatLift9x7
008 *
009 * Description:             An analyzing wavelet filter implementing the
010 *                          lifting 9x7 transform.
011 *
012 *
013 *
014 * COPYRIGHT:
015 *
016 * This software module was originally developed by Raphaël Grosbois and
017 * Diego Santa Cruz (Swiss Federal Institute of Technology-EPFL); Joel
018 * Askelöf (Ericsson Radio Systems AB); and Bertrand Berthelot, David
019 * Bouchard, Félix Henry, Gerard Mozelle and Patrice Onno (Canon Research
020 * Centre France S.A) in the course of development of the JPEG2000
021 * standard as specified by ISO/IEC 15444 (JPEG 2000 Standard). This
022 * software module is an implementation of a part of the JPEG 2000
023 * Standard. Swiss Federal Institute of Technology-EPFL, Ericsson Radio
024 * Systems AB and Canon Research Centre France S.A (collectively JJ2000
025 * Partners) agree not to assert against ISO/IEC and users of the JPEG
026 * 2000 Standard (Users) any of their rights under the copyright, not
027 * including other intellectual property rights, for this software module
028 * with respect to the usage by ISO/IEC and Users of this software module
029 * or modifications thereof for use in hardware or software products
030 * claiming conformance to the JPEG 2000 Standard. Those intending to use
031 * this software module in hardware or software products are advised that
032 * their use may infringe existing patents. The original developers of
033 * this software module, JJ2000 Partners and ISO/IEC assume no liability
034 * for use of this software module or modifications thereof. No license
035 * or right to this software module is granted for non JPEG 2000 Standard
036 * conforming products. JJ2000 Partners have full right to use this
037 * software module for his/her own purpose, assign or donate this
038 * software module to any third party and to inhibit third parties from
039 * using this software module for non JPEG 2000 Standard conforming
040 * products. This copyright notice must be included in all copies or
041 * derivative works of this software module.
042 *
043 * Copyright (c) 1999/2000 JJ2000 Partners.
044 * */
045package jj2000.j2k.wavelet.analysis;
046
047import jj2000.j2k.wavelet.*;
048import jj2000.j2k.image.*;
049import jj2000.j2k.*;
050import jj2000.j2k.codestream.writer.*;
051
052/**
053 * This class inherits from the analysis wavelet filter definition
054 * for int data. It implements the forward wavelet transform
055 * specifically for the 9x7 filter. The implementation is based on
056 * the lifting scheme.
057 *
058 * <P>See the AnWTFilter class for details such as
059 * normalization, how to split odd-length signals, etc. In particular,
060 * this method assumes that the low-pass coefficient is computed first.
061 *
062 * @see AnWTFilter
063 * @see AnWTFilterFloat
064 * */
065public class AnWTFilterFloatLift9x7 extends AnWTFilterFloat {
066
067    /** The low-pass synthesis filter of the 9x7 wavelet transform */
068    private final static float LPSynthesisFilter[] =
069    { -0.091272f, -0.057544f, 0.591272f, 1.115087f,
070      0.591272f, -0.057544f, -0.091272f};
071
072    /** The high-pass synthesis filter of the 9x7 wavelet transform */
073    private final static float HPSynthesisFilter[] =
074        { 0.026749f, 0.016864f, -0.078223f, -0.266864f,
075          0.602949f,
076          -0.266864f, -0.078223f, 0.016864f, 0.026749f };
077
078    /** The value of the first lifting step coefficient */
079    public final static float ALPHA = -1.586134342f;
080
081    /** The value of the second lifting step coefficient */
082    public final static float BETA = -0.05298011854f;
083
084    /** The value of the third lifting step coefficient */
085    public final static float GAMMA = 0.8829110762f;
086
087    /** The value of the fourth lifting step coefficient */
088    public final static float DELTA = 0.443568522f;
089
090    /** The value of the low-pass subband normalization factor */
091    public final static float KL = 0.8128930655f;//1.149604398f;
092
093    /** The value of the high-pass subband normalization factor */
094    public final static float KH = 1.230174106f;//0.8698644523f;
095
096    /**
097     * An implementation of the analyze_lpf() method that works on int
098     * data, for the forward 9x7 wavelet transform using the
099     * lifting scheme. See the general description of the analyze_lpf()
100     * method in the AnWTFilter class for more details.
101     *
102     * <P>The coefficients of the first lifting step are [ALPHA 1 ALPHA].
103     *
104     * <P>The coefficients of the second lifting step are [BETA 1 BETA].
105     *
106     * <P>The coefficients of the third lifting step are [GAMMA 1 GAMMA].
107     *
108     * <P>The coefficients of the fourth lifting step are [DELTA 1 DELTA].
109     *
110     * <P>The low-pass and high-pass subbands are normalized by respectively
111     * a factor of KL and a factor of KH
112     *
113     * @param inSig This is the array that contains the input
114     * signal.
115     *
116     * @param inOff This is the index in inSig of the first sample to
117     * filter.
118     *
119     * @param inLen This is the number of samples in the input signal
120     * to filter.
121     *
122     * @param inStep This is the step, or interleave factor, of the
123     * input signal samples in the inSig array.
124     *
125     * @param lowSig This is the array where the low-pass output
126     * signal is placed.
127     *
128     * @param lowOff This is the index in lowSig of the element where
129     * to put the first low-pass output sample.
130     *
131     * @param lowStep This is the step, or interleave factor, of the
132     * low-pass output samples in the lowSig array.
133     *
134     * @param highSig This is the array where the high-pass output
135     * signal is placed.
136     *
137     * @param highOff This is the index in highSig of the element where
138     * to put the first high-pass output sample.
139     *
140     * @param highStep This is the step, or interleave factor, of the
141     * high-pass output samples in the highSig array.
142     * */
143    public
144        void analyze_lpf(float inSig[], int inOff, int inLen, int inStep,
145                     float lowSig[], int lowOff, int lowStep,
146                     float highSig[], int highOff, int highStep) {
147        int i,maxi;
148        int iStep = 2 * inStep; //Subsampling in inSig
149        int ik;    //Indexing inSig
150        int lk;    //Indexing lowSig
151        int hk;    //Indexing highSig
152
153        // Generate intermediate high frequency subband
154
155        //Initialize counters
156        ik = inOff + inStep;
157        lk = lowOff;
158        hk = highOff;
159
160        //Apply first lifting step to each "inner" sample
161        for( i = 1, maxi = inLen-1; i < maxi; i += 2 ) {
162            highSig[hk] = inSig[ik] +
163                ALPHA*(inSig[ik-inStep] + inSig[ik+inStep]);
164
165            ik += iStep;
166            hk += highStep;
167        }
168
169        //Handle head boundary effect if input signal has even length
170        if(inLen % 2 == 0) {
171           highSig[hk] = inSig[ik] + 2*ALPHA*inSig[ik-inStep];
172        }
173
174        // Generate intermediate low frequency subband
175
176        //Initialize counters
177        ik = inOff;
178        lk = lowOff;
179        hk = highOff;
180
181        if(inLen>1) {
182            lowSig[lk] = inSig[ik] + 2*BETA*highSig[hk];
183        }
184        else {
185            lowSig[lk] = inSig[ik];
186        }
187
188        ik += iStep;
189        lk += lowStep;
190        hk += highStep;
191
192        //Apply lifting step to each "inner" sample
193        for( i = 2, maxi = inLen-1; i < maxi; i += 2 ) {
194            lowSig[lk] = inSig[ik] +
195                BETA*(highSig[hk-highStep] + highSig[hk]);
196
197            ik += iStep;
198            lk += lowStep;
199            hk += highStep;
200        }
201
202        //Handle head boundary effect if input signal has odd length
203        if((inLen % 2 == 1)&&(inLen>2)) {
204            lowSig[lk] =  inSig[ik] + 2*BETA*highSig[hk-highStep];
205        }
206
207        // Generate high frequency subband
208
209        //Initialize counters
210        lk = lowOff;
211        hk = highOff;
212
213        //Apply first lifting step to each "inner" sample
214        for(i = 1, maxi = inLen-1; i < maxi; i += 2)  {
215            highSig[hk] += GAMMA*(lowSig[lk] + lowSig[lk+lowStep]);
216
217            lk += lowStep;
218            hk += highStep;
219        }
220
221        //Handle head boundary effect if input signal has even length
222        if(inLen % 2 == 0) {
223            highSig[hk] += 2*GAMMA*lowSig[lk];
224        }
225
226        // Generate low frequency subband
227
228        //Initialize counters
229        lk = lowOff;
230        hk = highOff;
231
232        //Handle tail boundary effect
233        //If access the overlap then perform the lifting step
234        if(inLen>1){
235            lowSig[lk] += 2*DELTA*highSig[hk];
236        }
237
238        lk += lowStep;
239        hk += highStep;
240
241        //Apply lifting step to each "inner" sample
242        for(i = 2, maxi = inLen-1; i < maxi; i += 2) {
243            lowSig[lk] +=
244                DELTA*(highSig[hk - highStep] + highSig[hk]);
245
246            lk += lowStep;
247            hk += highStep;
248        }
249
250        //Handle head boundary effect if input signal has odd length
251        if((inLen % 2 == 1)&&(inLen>2)) {
252            lowSig[lk] +=  2*DELTA*highSig[hk-highStep];
253        }
254
255        // Normalize low and high frequency subbands
256
257        //Re-initialize counters
258        lk = lowOff;
259        hk = highOff;
260
261        //Normalize each sample
262        for( i=0 ; i<(inLen>>1); i++ ) {
263            lowSig[lk] *= KL;
264            highSig[hk] *= KH;
265            lk += lowStep;
266            hk += highStep;
267        }
268        //If the input signal has odd length then normalize the last low-pass
269        //coefficient (if input signal is length one filter is identity)
270        if( inLen%2==1 && inLen != 1) {
271            lowSig[lk] *= KL;
272        }
273    }
274
275    /**
276     * An implementation of the analyze_hpf() method that works on int
277     * data, for the forward 9x7 wavelet transform using the
278     * lifting scheme. See the general description of the analyze_hpf() method
279     * in the AnWTFilter class for more details.
280     *
281     * <P>The coefficients of the first lifting step are [ALPHA 1 ALPHA].
282     *
283     * <P>The coefficients of the second lifting step are [BETA 1 BETA].
284     *
285     * <P>The coefficients of the third lifting step are [GAMMA 1 GAMMA].
286     *
287     * <P>The coefficients of the fourth lifting step are [DELTA 1 DELTA].
288     *
289     * <P>The low-pass and high-pass subbands are normalized by respectively
290     * a factor of KL and a factor of KH
291     *
292     * @param inSig This is the array that contains the input
293     * signal.
294     *
295     * @param inOff This is the index in inSig of the first sample to
296     * filter.
297     *
298     * @param inLen This is the number of samples in the input signal
299     * to filter.
300     *
301     * @param inStep This is the step, or interleave factor, of the
302     * input signal samples in the inSig array.
303     *
304     * @param lowSig This is the array where the low-pass output
305     * signal is placed.
306     *
307     * @param lowOff This is the index in lowSig of the element where
308     * to put the first low-pass output sample.
309     *
310     * @param lowStep This is the step, or interleave factor, of the
311     * low-pass output samples in the lowSig array.
312     *
313     * @param highSig This is the array where the high-pass output
314     * signal is placed.
315     *
316     * @param highOff This is the index in highSig of the element where
317     * to put the first high-pass output sample.
318     *
319     * @param highStep This is the step, or interleave factor, of the
320     * high-pass output samples in the highSig array.
321     *
322     * @see AnWTFilter#analyze_hpf
323     * */
324    public void analyze_hpf(float inSig[], int inOff, int inLen, int inStep,
325                    float lowSig[], int lowOff, int lowStep,
326                    float highSig[], int highOff, int highStep) {
327
328        int i,maxi;
329        int iStep = 2 * inStep; //Subsampling in inSig
330        int ik;    //Indexing inSig
331        int lk;    //Indexing lowSig
332        int hk;    //Indexing highSig
333
334        // Generate intermediate high frequency subband
335
336        //Initialize counters
337        ik = inOff;
338        lk = lowOff;
339        hk = highOff;
340
341        if ( inLen>1 ) {
342            // apply symmetric extension.
343            highSig[hk] = inSig[ik] + 2*ALPHA*inSig[ik+inStep];
344        }
345        else {
346            // Normalize for Nyquist gain
347            highSig[hk] = inSig[ik]*2;
348        }
349
350        ik += iStep;
351        hk += highStep;
352
353        //Apply first lifting step to each "inner" sample
354        for( i = 2 ; i < inLen-1 ; i += 2 ) {
355            highSig[hk] = inSig[ik] +
356                ALPHA*(inSig[ik-inStep] + inSig[ik+inStep]);
357            ik += iStep;
358            hk += highStep;
359        }
360
361        //If input signal has odd length then we perform the lifting step
362        // i.e. apply a symmetric extension.
363        if( (inLen%2==1) && (inLen>1) ) {
364            highSig[hk] = inSig[ik] + 2*ALPHA*inSig[ik-inStep];
365        }
366
367        // Generate intermediate low frequency subband
368
369        //Initialize counters
370        //ik = inOff + inStep;
371        ik = inOff + inStep;
372        lk = lowOff;
373        hk = highOff;
374
375        //Apply lifting step to each "inner" sample
376        // we are at the component boundary
377        for(i = 1; i < inLen-1; i += 2) {
378            lowSig[lk] = inSig[ik] +
379                BETA*(highSig[hk] + highSig[hk+highStep]);
380
381            ik += iStep;
382            lk += lowStep;
383            hk += highStep;
384        }
385        if ( inLen>1 && inLen%2==0 ) {
386            // symetric extension
387            lowSig[lk] = inSig[ik]+2*BETA*highSig[hk];
388        }
389
390        // Generate high frequency subband
391
392        //Initialize counters
393        lk = lowOff;
394        hk = highOff;
395
396        if ( inLen>1 ) {
397            // symmetric extension.
398            highSig[hk] += GAMMA*2*lowSig[lk];
399        }
400        //lk += lowStep;
401        hk += highStep;
402
403        //Apply first lifting step to each "inner" sample
404        for(i = 2 ; i < inLen-1 ; i += 2)  {
405            highSig[hk] += GAMMA*(lowSig[lk] + lowSig[lk+lowStep]);
406            lk += lowStep;
407            hk += highStep;
408        }
409
410        //Handle head boundary effect
411        if ( inLen>1 && inLen%2==1 ) {
412            // symmetric extension.
413            highSig[hk] += GAMMA*2*lowSig[lk];
414        }
415
416        // Generate low frequency subband
417
418        //Initialize counters
419        lk = lowOff;
420        hk = highOff;
421
422        // we are at the component boundary
423        for(i = 1 ; i < inLen-1; i += 2) {
424            lowSig[lk] += DELTA*(highSig[hk] + highSig[hk+highStep]);
425            lk += lowStep;
426            hk += highStep;
427        }
428
429        if ( inLen>1 && inLen%2==0 ) {
430            lowSig[lk] += DELTA*2*highSig[hk];
431        }
432
433        // Normalize low and high frequency subbands
434
435        //Re-initialize counters
436        lk = lowOff;
437        hk = highOff;
438
439        //Normalize each sample
440        for( i=0 ; i<(inLen>>1); i++ ) {
441            lowSig[lk] *= KL;
442            highSig[hk] *= KH;
443            lk += lowStep;
444            hk += highStep;
445        }
446        //If the input signal has odd length then normalize the last high-pass
447        //coefficient (if input signal is length one filter is identity)
448        if( inLen%2==1 && inLen != 1) {
449            highSig[hk] *= KH;
450        }
451    }
452
453    /**
454     * Returns the negative support of the low-pass analysis
455     * filter. That is the number of taps of the filter in the
456     * negative direction.
457     *
458     * @return 2
459     * */
460    public int getAnLowNegSupport() {
461        return 4;
462    }
463
464    /**
465     * Returns the positive support of the low-pass analysis
466     * filter. That is the number of taps of the filter in the
467     * negative direction.
468     *
469     * @return The number of taps of the low-pass analysis filter in
470     * the positive direction
471     * */
472    public int getAnLowPosSupport() {
473        return 4;
474    }
475
476    /**
477     * Returns the negative support of the high-pass analysis
478     * filter. That is the number of taps of the filter in the
479     * negative direction.
480     *
481     * @return The number of taps of the high-pass analysis filter in
482     * the negative direction
483     * */
484    public int getAnHighNegSupport() {
485        return 3;
486    }
487
488    /**
489     * Returns the positive support of the high-pass analysis
490     * filter. That is the number of taps of the filter in the
491     * negative direction.
492     *
493     * @return The number of taps of the high-pass analysis filter in
494     * the positive direction
495     * */
496    public int getAnHighPosSupport() {
497        return 3;
498    }
499
500    /**
501     * Returns the negative support of the low-pass synthesis
502     * filter. That is the number of taps of the filter in the
503     * negative direction.
504     *
505     * <P>A MORE PRECISE DEFINITION IS NEEDED
506     *
507     * @return The number of taps of the low-pass synthesis filter in
508     * the negative direction
509     * */
510    public int getSynLowNegSupport() {
511        return 3;
512    }
513
514    /**
515     * Returns the positive support of the low-pass synthesis
516     * filter. That is the number of taps of the filter in the
517     * negative direction.
518     *
519     * <P>A MORE PRECISE DEFINITION IS NEEDED
520     *
521     * @return The number of taps of the low-pass synthesis filter in
522     * the positive direction
523     * */
524    public int getSynLowPosSupport() {
525        return 3;
526    }
527
528    /**
529     * Returns the negative support of the high-pass synthesis
530     * filter. That is the number of taps of the filter in the
531     * negative direction.
532     *
533     * <P>A MORE PRECISE DEFINITION IS NEEDED
534     *
535     * @return The number of taps of the high-pass synthesis filter in
536     * the negative direction
537     * */
538    public int getSynHighNegSupport() {
539        return 4;
540    }
541
542    /**
543     * Returns the positive support of the high-pass synthesis
544     * filter. That is the number of taps of the filter in the
545     * negative direction.
546     *
547     * <P>A MORE PRECISE DEFINITION IS NEEDED
548     *
549     * @return The number of taps of the high-pass synthesis filter in
550     * the positive direction
551     * */
552    public int getSynHighPosSupport() {
553        return 4;
554    }
555
556    /**
557     * Returns the time-reversed low-pass synthesis waveform of the
558     * filter, which is the low-pass filter. This is the time-reversed
559     * impulse response of the low-pass synthesis filter. It is used
560     * to calculate the L2-norm of the synthesis basis functions for a
561     * particular subband (also called energy weight).
562     *
563     * <P>The returned array may not be modified (i.e. a reference to
564     * the internal array may be returned by the implementation of
565     * this method).
566     *
567     * @return The time-reversed low-pass synthesis waveform of the
568     * filter.
569     * */
570    public float[] getLPSynthesisFilter() {
571        return LPSynthesisFilter;
572    }
573
574    /**
575     * Returns the time-reversed high-pass synthesis waveform of the
576     * filter, which is the high-pass filter. This is the
577     * time-reversed impulse response of the high-pass synthesis
578     * filter. It is used to calculate the L2-norm of the synthesis
579     * basis functions for a particular subband (also called energy
580     * weight).
581     *
582     * <P>The returned array may not be modified (i.e. a reference to
583     * the internal array may be returned by the implementation of
584     * this method).
585     *
586     * @return The time-reversed high-pass synthesis waveform of the
587     * filter.
588     * */
589    public float[] getHPSynthesisFilter() {
590        return HPSynthesisFilter;
591    }
592
593    /**
594     * Returns the implementation type of this filter, as defined in
595     * this class, such as WT_FILTER_INT_LIFT, WT_FILTER_FLOAT_LIFT,
596     * WT_FILTER_FLOAT_CONVOL.
597     *
598     * @return WT_FILTER_INT_LIFT.
599     * */
600    public int getImplType() {
601        return WT_FILTER_FLOAT_LIFT;
602    }
603
604    /**
605     * Returns the reversibility of the filter. A filter is considered
606     * reversible if it is suitable for lossless coding.
607     *
608     * @return true since the 9x7 is reversible, provided the appropriate
609     * rounding is performed.
610     * */
611    public boolean isReversible() {
612        return false;
613    }
614
615    /**
616     * Returns true if the wavelet filter computes or uses the
617     * same "inner" subband coefficient as the full frame wavelet transform,
618     * and false otherwise. In particular, for block based transforms with
619     * reduced overlap, this method should return false. The term "inner"
620     * indicates that this applies only with respect to the coefficient that
621     * are not affected by image boundaries processings such as symmetric
622     * extension, since there is not reference method for this.
623     *
624     * <P>The result depends on the length of the allowed overlap when
625     * compared to the overlap required by the wavelet filter. It also
626     * depends on how overlap processing is implemented in the wavelet
627     * filter.
628     *
629     * @param tailOvrlp This is the number of samples in the input
630     * signal before the first sample to filter that can be used for
631     * overlap.
632     *
633     * @param headOvrlp This is the number of samples in the input
634     * signal after the last sample to filter that can be used for
635     * overlap.
636     *
637     * @param inLen This is the lenght of the input signal to filter.The
638     * required number of samples in the input signal after the last sample
639     * depends on the length of the input signal.
640     *
641     * @return true if both overlaps are greater than 2, and correct
642     * processing is applied in the analyze() method.
643     * */
644    public boolean isSameAsFullWT(int tailOvrlp, int headOvrlp, int inLen) {
645
646        //If the input signal has even length.
647        if( inLen % 2 == 0) {
648            if( tailOvrlp >= 4 && headOvrlp >= 3 ) return true;
649            else return false;
650        }
651        //Else if the input signal has odd length.
652        else {
653            if( tailOvrlp >= 4 && headOvrlp >= 4 ) return true;
654            else return false;
655        }
656    }
657
658    /**
659     * Tests if the 'obj' object is the same filter as this one. Two filters
660     * are the same if the same filter code should be output for both filters
661     * by the encodeFilterCode() method.
662     *
663     * <P>Currently the implementation of this method only tests if 'obj' is
664     * also of the class AnWTFilterFloatLift9x7
665     *
666     * @param The object against which to test inequality.
667     * */
668    public boolean equals(Object obj) {
669        // To spped up test, first test for reference equality
670        return obj == this ||
671            obj instanceof AnWTFilterFloatLift9x7;
672    }
673
674    /**
675     * Returns the type of filter used according to the FilterTypes
676     * interface(W9x7).
677     *
678     * @see FilterTypes
679     *
680     * @return The filter type.
681     * */
682    public int getFilterType(){
683        return FilterTypes.W9X7;
684    }
685
686    /** Debugging method */
687    public String toString(){
688        return "w9x7";
689    }
690}