t1.c

This document is intended to serve as an introduction to Wavelet processing through a set of Matlab

/*
 * Copyright (c) 2001-2002, David Janssens
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "t1.h"
#include "j2k.h"
#include "mqc.h"
#include "int.h"
#include "mct.h"
#include "dwt.h"
#include "int.h"
#include "fix.h"
#include <stdio.h>
#include <math.h>

#define T1_MAXCBLKW 64
#define T1_MAXCBLKH 64

#define T1_SIG_NE 0x0001
#define T1_SIG_SE 0x0002
#define T1_SIG_SW 0x0004
#define T1_SIG_NW 0x0008
#define T1_SIG_N 0x0010
#define T1_SIG_E 0x0020
#define T1_SIG_S 0x0040
#define T1_SIG_W 0x0080
#define T1_SIG_OTH (T1_SIG_N|T1_SIG_NE|T1_SIG_E|T1_SIG_SE|T1_SIG_S|T1_SIG_SW|T1_SIG_W|T1_SIG_NW)
#define T1_SIG_PRIM (T1_SIG_N|T1_SIG_E|T1_SIG_S|T1_SIG_W)

#define T1_SGN_N 0x0100
#define T1_SGN_E 0x0200
#define T1_SGN_S 0x0400
#define T1_SGN_W 0x0800
#define T1_SGN (T1_SGN_N|T1_SGN_E|T1_SGN_S|T1_SGN_W)

#define T1_SIG 0x1000
#define T1_REFINE 0x2000
#define T1_VISIT 0x4000

#define T1_NUMCTXS_AGG 1
#define T1_NUMCTXS_ZC 9
#define T1_NUMCTXS_MAG 3
#define T1_NUMCTXS_SC 5
#define T1_NUMCTXS_UNI 1

#define T1_CTXNO_AGG 0
#define T1_CTXNO_ZC (T1_CTXNO_AGG+T1_NUMCTXS_AGG)
#define T1_CTXNO_MAG (T1_CTXNO_ZC+T1_NUMCTXS_ZC)
#define T1_CTXNO_SC (T1_CTXNO_MAG+T1_NUMCTXS_MAG)
#define T1_CTXNO_UNI (T1_CTXNO_SC+T1_NUMCTXS_SC)
#define T1_NUMCTXS (T1_CTXNO_UNI+T1_NUMCTXS_UNI)

#define T1_NMSEDEC_BITS 7
#define T1_NMSEDEC_FRACBITS (T1_NMSEDEC_BITS-1)

int t1_lut_ctxno_zc[1024];
int t1_lut_ctxno_sc[256];
int t1_lut_ctxno_mag[4096];
int t1_lut_spb[256];
int t1_lut_nmsedec_sig[1<<T1_NMSEDEC_BITS];
int t1_lut_nmsedec_sig0[1<<T1_NMSEDEC_BITS];
int t1_lut_nmsedec_ref[1<<T1_NMSEDEC_BITS];
int t1_lut_nmsedec_ref0[1<<T1_NMSEDEC_BITS];

int t1_data[T1_MAXCBLKH][T1_MAXCBLKH];
int t1_flags[T1_MAXCBLKH+2][T1_MAXCBLKH+2];

int t1_getctxno_zc(int f, int orient) {
    return t1_lut_ctxno_zc[(orient<<8)|(f&T1_SIG_OTH)];
}

int t1_getctxno_sc(int f) {
    return t1_lut_ctxno_sc[(f&(T1_SIG_PRIM|T1_SGN))>>4];
}

int t1_getctxno_mag(int f) {
    return t1_lut_ctxno_mag[(f&T1_SIG_OTH)|(((f&T1_REFINE)!=0)<<11)];
}

int t1_getspb(int f) {
    return t1_lut_spb[(f&(T1_SIG_PRIM|T1_SGN))>>4];
}

int t1_getnmsedec_sig(int x, int bitpos) {
    if (bitpos>T1_NMSEDEC_FRACBITS) {
        return t1_lut_nmsedec_sig[(x>>(bitpos-T1_NMSEDEC_FRACBITS))&((1<<T1_NMSEDEC_BITS)-1)];
    } else {
        return t1_lut_nmsedec_sig0[x&((1<<T1_NMSEDEC_BITS)-1)];
    }
}

int t1_getnmsedec_ref(int x, int bitpos) {
    if (bitpos>T1_NMSEDEC_FRACBITS) {
        return t1_lut_nmsedec_ref[(x>>(bitpos-T1_NMSEDEC_FRACBITS))&((1<<T1_NMSEDEC_BITS)-1)];
    } else {
        return t1_lut_nmsedec_ref0[x&((1<<T1_NMSEDEC_BITS)-1)];
    }
}

void t1_updateflags(int *fp, int s) {
    int *np=fp-(T1_MAXCBLKW+2);
    int *sp=fp+(T1_MAXCBLKW+2);
    np[-1]|=T1_SIG_SE;
    np[1]|=T1_SIG_SW;
    sp[-1]|=T1_SIG_NE;
    sp[1]|=T1_SIG_NW;
    *np|=T1_SIG_S;
    *sp|=T1_SIG_N;
    fp[-1]|=T1_SIG_E;
    fp[1]|=T1_SIG_W;
    if (s) {
        *np|=T1_SGN_S;
        *sp|=T1_SGN_N;
        fp[-1]|=T1_SGN_E;
        fp[1]|=T1_SGN_W;
    }
}

void t1_enc_sigpass_step(int *fp, int *dp, int orient, int bpno, int one, int *nmsedec) {
    int v;
    if ((*fp&T1_SIG_OTH) && !(*fp&(T1_SIG|T1_VISIT))) {
        mqc_setcurctx(t1_getctxno_zc(*fp, orient));
        v=int_abs(*dp)&one?1:0;
        mqc_encode(v);
        if (v) {
            *nmsedec+=t1_getnmsedec_sig(int_abs(*dp), bpno+T1_NMSEDEC_FRACBITS);
            mqc_setcurctx(t1_getctxno_sc(*fp));
            v=*dp<0?1:0;
            mqc_encode(v^t1_getspb(*fp));
            t1_updateflags(fp, v);
            *fp|=T1_SIG;
        }
        *fp|=T1_VISIT;
    }
}

void t1_dec_sigpass_step(int *fp, int *dp, int orient, int oneplushalf) {
    int v;
    if ((*fp&T1_SIG_OTH) && !(*fp&(T1_SIG|T1_VISIT))) {
        mqc_setcurctx(t1_getctxno_zc(*fp, orient));
        if (mqc_decode()) {
            mqc_setcurctx(t1_getctxno_sc(*fp));
            v=mqc_decode()^t1_getspb(*fp);
            *dp=v?-oneplushalf:oneplushalf;
            t1_updateflags(fp, v);
            *fp|=T1_SIG;
        }
        *fp|=T1_VISIT;
    }
}

void t1_enc_sigpass(int w, int h, int bpno, int orient, int *nmsedec) {
    int i, j, k, one;
    *nmsedec=0;
    one=1<<(bpno+T1_NMSEDEC_FRACBITS);
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            for (j=k; j<k+4 && j<h; j++) {
                t1_enc_sigpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], orient, bpno, one, nmsedec);
            }
        }
    }
}

void t1_dec_sigpass(int w, int h, int bpno, int orient) {
    int i, j, k, one, half, oneplushalf;
    one=1<<bpno;
    half=one>>1;
    oneplushalf=one|half;
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            for (j=k; j<k+4 && j<h; j++) {
                t1_dec_sigpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], orient, oneplushalf);
            }
        }
    }
}

void t1_enc_refpass_step(int *fp, int *dp, int bpno, int one, int *nmsedec) {
    int v;
    if ((*fp&(T1_SIG|T1_VISIT))==T1_SIG) {
        *nmsedec+=t1_getnmsedec_ref(int_abs(*dp), bpno+T1_NMSEDEC_FRACBITS);
        mqc_setcurctx(t1_getctxno_mag(*fp));
        v=int_abs(*dp)&one?1:0;
        mqc_encode(v);
        *fp|=T1_REFINE;
    }
}

void t1_dec_refpass_step(int *fp, int *dp, int poshalf, int neghalf) {
    int v, t;
    if ((*fp&(T1_SIG|T1_VISIT))==T1_SIG) {
        mqc_setcurctx(t1_getctxno_mag(*fp));
        v=mqc_decode();
        t=v?poshalf:neghalf;
        *dp+=*dp<0?-t:t;
        *fp|=T1_REFINE;
    }
}

void t1_enc_refpass(int w, int h, int bpno, int *nmsedec) {
    int i, j, k, one;
    *nmsedec=0;
    one=1<<(bpno+T1_NMSEDEC_FRACBITS);
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            for (j=k; j<k+4 && j<h; j++) {
                t1_enc_refpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], bpno, one, nmsedec);
            }
        }
    }
}

void t1_dec_refpass(int w, int h, int bpno) {
    int i, j, k, one, poshalf, neghalf;
    one=1<<bpno;
    poshalf=one>>1;
    neghalf=bpno>0?-poshalf:-1;
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            for (j=k; j<k+4 && j<h; j++) {
                t1_dec_refpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], poshalf, neghalf);
            }
        }
    }
}

void t1_enc_clnpass_step(int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, int partial) {
    int v;
    if (partial) goto label_partial;
    if (!(*fp&(T1_SIG|T1_VISIT))) {
        mqc_setcurctx(t1_getctxno_zc(*fp, orient));
        v=int_abs(*dp)&one?1:0;
        mqc_encode(v);
        if (v) {
label_partial:
            *nmsedec+=t1_getnmsedec_sig(int_abs(*dp), bpno+T1_NMSEDEC_FRACBITS);
            mqc_setcurctx(t1_getctxno_sc(*fp));
            v=*dp<0?1:0;
            mqc_encode(v^t1_getspb(*fp));
            t1_updateflags(fp, v);
            *fp|=T1_SIG;
        }
    }
    *fp&=~T1_VISIT;
}

void t1_dec_clnpass_step(int *fp, int *dp, int orient, int oneplushalf, int partial) {
    int v;
    if (partial) goto label_partial;
    if (!(*fp&(T1_SIG|T1_VISIT))) {
        mqc_setcurctx(t1_getctxno_zc(*fp, orient));
        if (mqc_decode()) {
label_partial:
            mqc_setcurctx(t1_getctxno_sc(*fp));
            v=mqc_decode()^t1_getspb(*fp);
            *dp=v?-oneplushalf:oneplushalf;
            t1_updateflags(fp, v);
            *fp|=T1_SIG;
        }
    }
    *fp&=~T1_VISIT;
}

void t1_enc_clnpass(int w, int h, int bpno, int orient, int *nmsedec) {
    int i, j, k, one, agg, runlen;
    *nmsedec=0;
    one=1<<(bpno+T1_NMSEDEC_FRACBITS);
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            if (k+3<h) {
                agg=!(t1_flags[1+k][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+1][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+2][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+3][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH));
            } else {
                agg=0;
            }
            if (agg) {
                for (runlen=0; runlen<4; runlen++) {
                    if (int_abs(t1_data[k+runlen][i])&one) break;
                }
                mqc_setcurctx(T1_CTXNO_AGG);
                mqc_encode(runlen!=4);
                if (runlen==4) {
                    continue;
                }
                mqc_setcurctx(T1_CTXNO_UNI);
                mqc_encode(runlen>>1);
                mqc_encode(runlen&1);
            } else {
                runlen=0;
            }
            for (j=k+runlen; j<k+4 && j<h; j++) {
                t1_enc_clnpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], orient, bpno, one, nmsedec, agg && (j==k+runlen));
            }
        }
    }
}

void t1_dec_clnpass(int w, int h, int bpno, int orient, int segsym) {
    int i, j, k, one, half, oneplushalf, agg, runlen;
    one=1<<bpno;
    half=one>>1;
    oneplushalf=one|half;
    for (k=0; k<h; k+=4) {
        for (i=0; i<w; i++) {
            if (k+3<h) {
                agg=!(t1_flags[1+k][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+1][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+2][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH) ||
                        t1_flags[1+k+3][1+i]&(T1_SIG|T1_VISIT|T1_SIG_OTH));
            } else {
                agg=0;
            }
            if (agg) {
                mqc_setcurctx(T1_CTXNO_AGG);
                if (!mqc_decode()) {
                    continue;
                }
                mqc_setcurctx(T1_CTXNO_UNI);
                runlen=mqc_decode();
                runlen=(runlen<<1)|mqc_decode();
            } else {
                runlen=0;
            }
            for (j=k+runlen; j<k+4 && j<h; j++) {
                t1_dec_clnpass_step(&t1_flags[1+j][1+i], &t1_data[j][i], orient, oneplushalf, agg && (j==k+runlen));