t1.c

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

            }
        }
    }
    if (segsym) {
        int v=0;
        mqc_setcurctx(T1_CTXNO_UNI);
        v=mqc_decode();
        v=(v<<1)|mqc_decode();
        v=(v<<1)|mqc_decode();
        v=(v<<1)|mqc_decode();
        if (v!=0xa) {
      /*             fprintf(stderr, "warning: bad segmentation symbol\n"); */
        }
    }
}

double t1_getwmsedec(int nmsedec, int compno, int level, int orient, int bpno, int qmfbid, double stepsize) {
    double w1, w2, wmsedec;
    if (qmfbid==1) {
        w1=mct_getnorm(compno);
        w2=dwt_getnorm(level, orient);
    } else if (qmfbid==0) {
        w1=mct_getnorm_real(compno);
        w2=dwt_getnorm_real(level, orient);
    }
    wmsedec=w1*w2*stepsize*(1<<bpno);
    wmsedec*=wmsedec*nmsedec/8192.0;
    return wmsedec;
}

void t1_encode_cblk(tcd_cblk_t *cblk, int orient, int compno, int level, int qmfbid, double stepsize) {
    int i, j;
    int w, h;
    int passno;
    int bpno, passtype;
    int max;
    int nmsedec;
    double cumwmsedec=0;

    w=cblk->x1-cblk->x0;
    h=cblk->y1-cblk->y0;

    max=0;
    for (j=0; j<h; j++) {
        for (i=0; i<w; i++) {
            max=int_max(max, int_abs(t1_data[j][i]));
        }
    }
    cblk->numbps=max?(int_floorlog2(max)+1)-T1_NMSEDEC_FRACBITS:0;

    for (i=0; i<sizeof(t1_flags)/sizeof(int); i++) ((int*)t1_flags)[i]=0;

    bpno=cblk->numbps-1;
    passtype=2;

    mqc_resetstates();
    mqc_setstate(T1_CTXNO_UNI, 0, 46);
    mqc_setstate(T1_CTXNO_AGG, 0, 3);
    mqc_setstate(T1_CTXNO_ZC, 0, 4);
    mqc_init_enc(cblk->data);

    for (passno=0; bpno>=0; passno++) {
        tcd_pass_t *pass=&cblk->passes[passno];
        switch (passtype) {
            case 0:
                t1_enc_sigpass(w, h, bpno, orient, &nmsedec);
                break;
            case 1:
                t1_enc_refpass(w, h, bpno, &nmsedec);
                break;
            case 2:
                t1_enc_clnpass(w, h, bpno, orient, &nmsedec);
                break;
        }
        cumwmsedec+=t1_getwmsedec(nmsedec, compno, level, orient, bpno, qmfbid, stepsize);
        pass->distortiondec=cumwmsedec;
        pass->rate=mqc_numbytes()+3; /*  FIXME */
        if (++passtype==3) {
            passtype=0;
            bpno--;
        }
    }
    mqc_flush();

    cblk->totalpasses=passno;
}

void t1_decode_cblk(tcd_cblk_t *cblk, int orient, int roishift, int cblksty) {
    int i;
    int w, h;
    int bpno, passtype;
    int segno, passno;

    for (i=0; i<sizeof(t1_data)/sizeof(int); i++) ((int*)t1_data)[i]=0;
    for (i=0; i<sizeof(t1_flags)/sizeof(int); i++) ((int*)t1_flags)[i]=0;

    w=cblk->x1-cblk->x0;
    h=cblk->y1-cblk->y0;
    bpno=roishift+cblk->numbps-1;
    passtype=2;

    mqc_resetstates();
    mqc_setstate(T1_CTXNO_UNI, 0, 46);
    mqc_setstate(T1_CTXNO_AGG, 0, 3);
    mqc_setstate(T1_CTXNO_ZC, 0, 4);

    for (segno=0; segno<cblk->numsegs; segno++) {
        tcd_seg_t *seg=&cblk->segs[segno];
        /* fprintf(stderr, "t1_decode_cblk: segno=%d, seg->len=%d, seg->numpasses=%d\n", segno, seg->len, seg->numpasses); */
        mqc_init_dec(seg->data, seg->len);
        for (passno=0; passno<seg->numpasses; passno++) {
            switch (passtype) {
                case 0:
                    t1_dec_sigpass(w, h, bpno, orient);
                    break;
                case 1:
                    t1_dec_refpass(w, h, bpno);
                    break;
                case 2:
                    t1_dec_clnpass(w, h, bpno, orient, cblksty&J2K_CCP_CBLKSTY_SEGSYM);
                    break;
            }
            if (++passtype==3) {
                passtype=0;
                bpno--;
            }
        }
    }
}

void t1_encode_cblks(tcd_tile_t *tile, j2k_tcp_t *tcp) {
    int compno, resno, bandno, precno, cblkno;
    for (compno=0; compno<tile->numcomps; compno++) {
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        for (resno=0; resno<tilec->numresolutions; resno++) {
            tcd_resolution_t *res=&tilec->resolutions[resno];
            for (bandno=0; bandno<res->numbands; bandno++) {
                tcd_band_t *band=&res->bands[bandno];
                for (precno=0; precno<res->pw*res->ph; precno++) {
                    tcd_precinct_t *prc=&band->precincts[precno];
                    for (cblkno=0; cblkno<prc->cw*prc->ch; cblkno++) {
                        int x, y, i, j, orient;
                        tcd_cblk_t *cblk=&prc->cblks[cblkno];
                        if (band->bandno==0) {
                            x=res->x0+cblk->x0;
                            y=res->y0+cblk->y0;
                        } else if (band->bandno==1) {
                            x=res->x1-(band->x1-band->x0)+cblk->x0;
                            y=res->y0+cblk->y0;
                        } else if (band->bandno==2) {
                            x=res->x0+cblk->x0;
                            y=res->y1-(band->y1-band->y0)+cblk->y0;
                        } else if (band->bandno==3) {
                            x=res->x1-(band->x1-band->x0)+cblk->x0;
                            y=res->y1-(band->y1-band->y0)+cblk->y0;
                        }
                        if (tcp->tccps[compno].qmfbid==1) {
                            for (j=0; j<cblk->y1-cblk->y0; j++) {
                                for (i=0; i<cblk->x1-cblk->x0; i++) {
                                    t1_data[j][i]=tilec->data[x+i+(y+j)*(tilec->x1-tilec->x0)]<<T1_NMSEDEC_FRACBITS;
                                }
                            }
                        } else if (tcp->tccps[compno].qmfbid==0) {
                            for (j=0; j<cblk->y1-cblk->y0; j++) {
                                for (i=0; i<cblk->x1-cblk->x0; i++) {
                                    t1_data[j][i]=fix_mul(tilec->data[x+i+(y+j)*(tilec->x1-tilec->x0)], 8192*8192/band->stepsize)>>(13-T1_NMSEDEC_FRACBITS);
                                }
                            }
                        }
                        orient=band->bandno; /*  FIXME */
                        if (orient==2) orient=1;
                        else if (orient==1) orient=2;
                        t1_encode_cblk(cblk, orient, compno, tilec->numresolutions-1-resno, tcp->tccps[compno].qmfbid, band->stepsize);
                    }
                }
            }
        }
    }
}

void t1_decode_cblks(tcd_tile_t *tile, j2k_tcp_t *tcp) {
    int compno, resno, bandno, precno, cblkno;
    for (compno=0; compno<tile->numcomps; compno++) {
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        for (resno=0; resno<tilec->numresolutions; resno++) {
            tcd_resolution_t *res=&tilec->resolutions[resno];
            for (bandno=0; bandno<res->numbands; bandno++) {
                tcd_band_t *band=&res->bands[bandno];
                for (precno=0; precno<res->pw*res->ph; precno++) {
                    tcd_precinct_t *prc=&band->precincts[precno];
                    for (cblkno=0; cblkno<prc->cw*prc->ch; cblkno++) {
                        int x, y, i, j, orient;
                        tcd_cblk_t *cblk=&prc->cblks[cblkno];
                        orient=band->bandno; /*  FIXME */
                        if (orient==2) orient=1;
                        else if (orient==1) orient=2;
                        t1_decode_cblk(cblk, orient, tcp->tccps[compno].roishift, tcp->tccps[compno].cblksty);
                        if (band->bandno==0) {
                            x=cblk->x0-band->x0;
                            y=cblk->y0-band->y0;
                        } else if (band->bandno==1) {
                            tcd_resolution_t *pres=&tilec->resolutions[resno-1];
                            x=pres->x1-pres->x0+cblk->x0-band->x0;
                            y=cblk->y0-band->y0;
                        } else if (band->bandno==2) {
                            tcd_resolution_t *pres=&tilec->resolutions[resno-1];
                            x=cblk->x0-band->x0;
                            y=pres->y1-pres->y0+cblk->y0-band->y0;
                        } else if (band->bandno==3) {
                            tcd_resolution_t *pres=&tilec->resolutions[resno-1];
                            x=pres->x1-pres->x0+cblk->x0-band->x0;
                            y=pres->y1-pres->y0+cblk->y0-band->y0;
                        }
                        if (tcp->tccps[compno].roishift) {
                            int thresh, val, mag;
                            thresh=1<<tcp->tccps[compno].roishift;
                            for (j=0; j<cblk->y1-cblk->y0; j++) {
                                for (i=0; i<cblk->x1-cblk->x0; i++) {
                                    val=t1_data[j][i];
                                    mag=int_abs(val);
                                    if (mag>=thresh) {
                                        mag>>=tcp->tccps[compno].roishift;
                                        t1_data[j][i]=val<0?-mag:mag;
                                    }
                                }
                            }
                        }
                        if (tcp->tccps[compno].qmfbid==1) {
                            for (j=0; j<cblk->y1-cblk->y0; j++) {
                                for (i=0; i<cblk->x1-cblk->x0; i++) {
                                    tilec->data[x+i+(y+j)*(tilec->x1-tilec->x0)]=t1_data[j][i];
                                }
                            }
                        } else if (tcp->tccps[compno].qmfbid==0) {
                            for (j=0; j<cblk->y1-cblk->y0; j++) {
                                for (i=0; i<cblk->x1-cblk->x0; i++) {
                                    if (t1_data[j][i]==0) {
                                        tilec->data[x+i+(y+j)*(tilec->x1-tilec->x0)]=0;
                                    } else {
                                        tilec->data[x+i+(y+j)*(tilec->x1-tilec->x0)]=fix_mul(t1_data[j][i]<<13, band->stepsize);
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

int t1_init_ctxno_zc(int f, int orient) {
    int h,v,d,n,t,hv;
    n=0;
    h=((f&T1_SIG_W)!=0)+((f&T1_SIG_E)!=0);
    v=((f&T1_SIG_N)!=0)+((f&T1_SIG_S)!=0);
    d=((f&T1_SIG_NW)!=0)+((f&T1_SIG_NE)!=0)+((f&T1_SIG_SE)!=0)+((f&T1_SIG_SW)!=0);
    switch(orient){
    case 2:
      t=h;h=v;v=t;
    case 0:
    case 1:
      if(!h){
    if(!v){
      if(!d){
        n=0;}
      else if(d==1){
        n=1;}
      else{
        n=2;}
    }
    else if(v==1){
      n=3;
    }
    else {
      n=4;}
      }
      else if(h==1){
    if(!v){
      if(!d){
        n=5;}
      else{
        n=6;}
    }
    else{
      n=7;}
      }
      else{
    n=8;}
      break;
    case 3:
      hv=h+v;
      if(!d){
    if(!hv){
      n=0;}
    else if(hv==1){
      n=1;}
    else{
      n=2;}
      }
      else if(d==1){
    if(!hv){n=3;}
    else if(hv==1){
      n=4;}
    else{n=5;}
      }
      else if(d==2)
    {if(!hv)
      {n=6;}
    else{n=7;}
    }
      else{n=8;}
      break;}
    return T1_CTXNO_ZC+n;
}

int t1_init_ctxno_sc(int f) {
    int hc,vc,n;n=0;hc=int_min(((f&(T1_SIG_E|T1_SGN_E))==T1_SIG_E)+((f&(T1_SIG_W|T1_SGN_W))==T1_SIG_W),1)-int_min(((f&(T1_SIG_E|T1_SGN_E))==(T1_SIG_E|T1_SGN_E))+((f&(T1_SIG_W|T1_SGN_W))==(T1_SIG_W|T1_SGN_W)),1);vc=int_min(((f&(T1_SIG_N|T1_SGN_N))==T1_SIG_N)+((f&(T1_SIG_S|T1_SGN_S))==T1_SIG_S),1)-int_min(((f&(T1_SIG_N|T1_SGN_N))==(T1_SIG_N|T1_SGN_N))+((f&(T1_SIG_S|T1_SGN_S))==(T1_SIG_S|T1_SGN_S)),1);if(hc<0){hc=-hc;vc=-vc;}if(!hc){if(vc==-1){n=1;}else if(!vc){n=0;}else{n=1;}}else if(hc==1){if(vc==-1){n=2;}else if(!vc){n=3;}else{n=4;}}return T1_CTXNO_SC+n;
}

int t1_init_ctxno_mag(int f) {
    int n;if(!(f&T1_REFINE)){n=(f&(T1_SIG_OTH))?1:0;}else{n=2;}return T1_CTXNO_MAG+n;
}

int t1_init_spb(int f) {
    int hc,vc,n;hc=int_min(((f&(T1_SIG_E|T1_SGN_E))==T1_SIG_E)+((f&(T1_SIG_W|T1_SGN_W))==T1_SIG_W),1)-int_min(((f&(T1_SIG_E|T1_SGN_E))==(T1_SIG_E|T1_SGN_E))+((f&(T1_SIG_W|T1_SGN_W))==(T1_SIG_W|T1_SGN_W)),1);vc=int_min(((f&(T1_SIG_N|T1_SGN_N))==T1_SIG_N)+((f&(T1_SIG_S|T1_SGN_S))==T1_SIG_S),1)-int_min(((f&(T1_SIG_N|T1_SGN_N))==(T1_SIG_N|T1_SGN_N))+((f&(T1_SIG_S|T1_SGN_S))==(T1_SIG_S|T1_SGN_S)),1);if(!hc&&!vc){n=0;}else{n=(!(hc>0||(!hc&&vc>0)));}return n;
}

void t1_init_luts() {
    int i,j;double u,v,t;for(j=0;j<4;j++){for(i=0;i<256;++i){t1_lut_ctxno_zc[(j<<8)|i]=t1_init_ctxno_zc(i,j);}}for(i=0;i<256;i++){t1_lut_ctxno_sc[i]=t1_init_ctxno_sc(i<<4);}for(j=0;j<2;j++){for(i=0;i<2048;++i){t1_lut_ctxno_mag[(j<<11)+i]=t1_init_ctxno_mag((j?T1_REFINE:0)|i);}}for(i=0;i<256;++i){t1_lut_spb[i]=t1_init_spb(i<<4);}
    /*  FIXME FIXME FIXME */
    /* printf("nmsedec luts:\n"); */
    for (i=0; i<(1<<T1_NMSEDEC_BITS); i++) {
        t=i/pow(2,T1_NMSEDEC_FRACBITS);
        u=t;
        v=t-1.5;
        t1_lut_nmsedec_sig[i]=int_max(0,(int)(floor((u*u-v*v)*pow(2,T1_NMSEDEC_FRACBITS)+ 0.5)/pow(2,T1_NMSEDEC_FRACBITS)*8192.0));
        t1_lut_nmsedec_sig0[i]=int_max(0,(int)(floor((u*u)*pow(2,T1_NMSEDEC_FRACBITS)+0.5)/pow(2,T1_NMSEDEC_FRACBITS)*8192.0));
        u=t-1.0;
        if (i&(1<<(T1_NMSEDEC_BITS-1))) {
            v=t-1.5;
        } else {
            v=t-0.5;
        }
        t1_lut_nmsedec_ref[i]=int_max(0,(int)(floor((u*u-v*v)*pow(2,T1_NMSEDEC_FRACBITS)+0.5)/pow(2,T1_NMSEDEC_FRACBITS)*8192.0));
        t1_lut_nmsedec_ref0[i]=int_max(0,(int)(floor((u*u)*pow(2,T1_NMSEDEC_FRACBITS)+0.5)/pow(2,T1_NMSEDEC_FRACBITS)*8192.0));
        /* printf("%d: %d %d %d %d\n", i, t1_lut_nmsedec_sig[i], t1_lut_nmsedec_sig0[i], t1_lut_nmsedec_ref[i], t1_lut_nmsedec_ref0[i]); */
    }
}