tcd.c

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

    max=0;
    for (compno=0; compno<tcd_tile->numcomps; compno++) {
        tcd_tilecomp_t *tilec=&tcd_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++) {
                        tcd_cblk_t *cblk=&prc->cblks[cblkno];
                        for (passno=0; passno<cblk->totalpasses; passno++) {
                            tcd_pass_t *pass=&cblk->passes[passno];
                            int dr;
                            double dd, rdslope;
                            if (passno==0) {
                                dr=pass->rate;
                                dd=pass->distortiondec;
                            } else {
                                dr=pass->rate-cblk->passes[passno-1].rate;
                                dd=pass->distortiondec-cblk->passes[passno-1].distortiondec;
                            }
                            if (dr==0) {
                                continue;
                            }
                            rdslope=dd/dr;
                            if (rdslope<min) {
                                min=rdslope;
                            }
                            if (rdslope>max) {
                                max=rdslope;
                            }
                        }
                    }
                }
            }
        }
    }
    for (layno=0; layno<tcd_tcp->numlayers; layno++) {
        volatile double lo=min;
        volatile double hi=max;
        volatile int success=0;
        volatile int maxlen=int_min(tcd_tcp->rates[layno], len);
        volatile double goodthresh;
        volatile int goodlen;
        volatile jmp_buf oldenv;
        volatile int i;
        memcpy((void*)oldenv, j2k_error, sizeof(jmp_buf));
        for (i=0; i<32; i++) {
            volatile double thresh=(lo+hi)/2;
            int l;
            tcd_makelayer(layno, thresh, 0);
            if (setjmp(j2k_error)) {
                lo=thresh;
                continue;
            }
            l=t2_encode_packets(tcd_img, tcd_cp, tcd_tileno, tcd_tile, layno+1, dest, maxlen);
        /*             fprintf(stderr, "rate alloc: len=%d, max=%d\n", l, maxlen); */
            hi=thresh;
            success=1;
            goodthresh=thresh;
            goodlen=l;
        }
        memcpy(j2k_error, (void*)oldenv, sizeof(jmp_buf));
        if (!success) {
            longjmp(j2k_error, 1);
        }
        tcd_makelayer(layno, goodthresh, 1);
    }
}

int tcd_encode_tile(int tileno, unsigned char *dest, int len) {


    int taille_x = (int) (*( (int *)taille_x_e));
    int taille_y = (int) (*( (int *)taille_y_e));


    int compno;
    int l;
#ifndef TIME_H_NOT_AVAILABLE
    clock_t time1, time2, time3, time4, time5, time6, time7;
#endif
    tcd_tile_t *tile;
    tcd_tileno=tileno;
    tcd_tile=&tcd_image.tiles[tileno];
    tcd_tcp=&tcd_cp->tcps[tileno];
    tile=tcd_tile;

#ifndef TIME_H_NOT_AVAILABLE
    time7=clock();

    time1=clock();
#endif
    for (compno=0; compno<tile->numcomps; compno++) {
        int i, j;
        int tw, w;
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        /*
          if signed, no adjustment is performed
          if unsigned, subtracting the median value : for example, 2 bits [0..3] -> [-2..1]
        */
        int adjust=tcd_img->comps[compno].sgnd?0:1<<(tcd_img->comps[compno].prec-1);
        tw=tilec->x1-tilec->x0;
        w=int_ceildiv(tcd_img->x1-tcd_img->x0, tcd_img->comps[compno].dx);
        for (j=tilec->y0; j<tilec->y1; j++) {
            for (i=tilec->x0; i<tilec->x1; i++) {
	        #ifndef ENCODE_ONLY
                if (tcd_tcp->tccps[compno].qmfbid==1) {
                    tilec->data[i-tilec->x0+(j-tilec->y0)*tw]=tcd_img->comps[compno].data[i+j*w]-adjust;
                } else if (tcd_tcp->tccps[compno].qmfbid==0) {
                    /* conversion to FIXED is performed for real calculations only */
                    tilec->data[i-tilec->x0+(j-tilec->y0)*tw]=(tcd_img->comps[compno].data[i+j*w]-adjust)<<13;
		}
		#endif
		#ifdef ENCODE_ONLY
			tilec->data[i-tilec->x0+(j-tilec->y0)*tw]= (int) ((*((double*)matlab_in_data + i + j*taille_x))*(1<<13));
		#endif
	   }
        }
    }
#ifndef TIME_H_NOT_AVAILABLE
    time1=clock()-time1;

    time2=clock();
#endif
#ifndef ENCODE_ONLY
    if (tcd_tcp->mct) {
        if (tcd_tcp->tccps[0].qmfbid==0) {
            mct_encode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, (tile->comps[0].x1-tile->comps[0].x0)*(tile->comps[0].y1-tile->comps[0].y0));
        } else {
            mct_encode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, (tile->comps[0].x1-tile->comps[0].x0)*(tile->comps[0].y1-tile->comps[0].y0));
        }
    }
#endif
#ifndef TIME_H_NOT_AVAILABLE
    time2=clock()-time2;

    time3=clock();
#endif
#ifndef ENCODE_ONLY
    for (compno=0; compno<tile->numcomps; compno++) { 
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        if (tcd_tcp->tccps[compno].qmfbid==1) {
            dwt_encode(tilec->data, tilec->x1-tilec->x0, tilec->y1-tilec->y0, tilec->numresolutions-1);
        } else if (tcd_tcp->tccps[compno].qmfbid==0)
        {	    dwt_encode_real(tilec->data, tilec->x1-tilec->x0, tilec->y1-tilec->y0, tilec->numresolutions-1);
        }
		    memcpy((int*) matlab_ex_data, (int*) tilec->data ,taille_x*taille_y*sizeof(int));
    }
#endif
#ifndef TIME_H_NOT_AVAILABLE
    time3=clock()-time3;

    time4=clock();
#endif
    t1_init_luts();
/*	fprintf(stderr, "taille des codeblocks : %d %d, type codage codeblocks : %d\n", tcd_tcp->tccps->cblkw, tcd_tcp->tccps->cblkh, tcd_tcp->tccps->cblksty);
	fprintf(stderr, "taille des precincts : %d %d\n", tcd_tcp->tccps->prcw , tcd_tcp->tccps->prch);
	fprintf(stderr, "Pas de quantif.(exposant/mantisse) : %d %d\n", tcd_tcp->tccps->stepsizes->expn,tcd_tcp->tccps->stepsizes->mant);
*/
    t1_encode_cblks(tile, tcd_tcp);
#ifndef TIME_H_NOT_AVAILABLE
    time4=clock()-time4;

    time5=clock();
#endif
    tcd_rateallocate(dest, len);
#ifndef TIME_H_NOT_AVAILABLE
    time5=clock()-time5;

    time6=clock();
#endif
    l=t2_encode_packets(tcd_img, tcd_cp, tileno, tile, tcd_tcp->numlayers, dest, len);
#ifndef TIME_H_NOT_AVAILABLE
    time6=clock()-time6;

    time7=clock()-time7;
#endif

    /*     printf("tile encoding times:\n"); */
    /*     printf("img->tile: %d.%.3d s\n", time1/CLOCKS_PER_SEC, (time1%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*       printf("mct:       %d.%.3d s\n", time2/CLOCKS_PER_SEC, (time2%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("dwt:       %d.%.3d s\n", time3/CLOCKS_PER_SEC, (time3%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("tier 1:    %d.%.3d s\n", time4/CLOCKS_PER_SEC, (time4%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("ratealloc: %d.%.3d s\n", time5/CLOCKS_PER_SEC, (time5%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("tier 2:    %d.%.3d s\n", time6/CLOCKS_PER_SEC, (time6%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("total:     %d.%.3d s\n", time7/CLOCKS_PER_SEC, (time7%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */

    return l;
}

int tcd_decode_tile(unsigned char *src, int len, int tileno) {

    #ifdef ENCODE_ONLY
    int taille_x = (int) (*( (int *)taille_x_e));
    int taille_y = (int) (*( (int *)taille_y_e));
    #endif

    int l;
    int compno;
    int eof=0;
    jmp_buf oldenv;
#ifndef TIME_H_NOT_AVAILABLE
    clock_t time1, time2, time3, time4, time5, time6;
#endif

    tcd_tile_t *tile;
    tcd_tileno=tileno;
    tcd_tile=&tcd_image.tiles[tileno];
    tcd_tcp=&tcd_cp->tcps[tileno];
    tile=tcd_tile;

#ifndef TIME_H_NOT_AVAILABLE
    time6=clock();

    time1=clock();
#endif
    memcpy(oldenv, j2k_error, sizeof(jmp_buf));
    if (setjmp(j2k_error)) {
        eof=1;
    /*         fprintf(stderr, "tcd_decode: incomplete bistream\n"); */
    } else {
        l=t2_decode_packets(src, len, tcd_img, tcd_cp, tileno, tile);
    }
    memcpy(j2k_error, oldenv, sizeof(jmp_buf));
#ifndef TIME_H_NOT_AVAILABLE
    time1=clock()-time1;

    time2=clock();
#endif
    t1_init_luts();
    t1_decode_cblks(tile, tcd_tcp);
#ifndef TIME_H_NOT_AVAILABLE
    time2=clock()-time2;

    time3=clock();
#endif
#ifndef ENCODE_ONLY
    for (compno=0; compno<tile->numcomps; compno++) {
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        if (tcd_tcp->tccps[compno].qmfbid==1) {
            dwt_decode(tilec->data, tilec->x1-tilec->x0, tilec->y1-tilec->y0, tilec->numresolutions-1);
        } else if (tcd_tcp->tccps[compno].qmfbid==0) {
            dwt_decode_real(tilec->data, tilec->x1-tilec->x0, tilec->y1-tilec->y0, tilec->numresolutions-1);
        }
    }
#endif
#ifndef TIME_H_NOT_AVAILABLE
    time3=clock()-time3;

    time4=clock();
#endif
#ifndef ENCODE_ONLY
    if (tcd_tcp->mct) {
        if (tcd_tcp->tccps[0].qmfbid==0) {
            mct_decode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, (tile->comps[0].x1-tile->comps[0].x0)*(tile->comps[0].y1-tile->comps[0].y0));
        } else {
            mct_decode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, (tile->comps[0].x1-tile->comps[0].x0)*(tile->comps[0].y1-tile->comps[0].y0));
        }
    }
#endif
#ifndef TIME_H_NOT_AVAILABLE
    time4=clock()-time4;

    time5=clock();
#endif

    for (compno=0; compno<tile->numcomps; compno++) {
        tcd_tilecomp_t *tilec=&tile->comps[compno];
        int adjust=tcd_img->comps[compno].sgnd?0:1<<(tcd_img->comps[compno].prec-1);
        int min=tcd_img->comps[compno].sgnd?-(1<<(tcd_img->comps[compno].prec-1)):0;
        int max=tcd_img->comps[compno].sgnd?(1<<(tcd_img->comps[compno].prec-1))-1:(1<<tcd_img->comps[compno].prec)-1;
        int tw=tilec->x1-tilec->x0;
        int w=int_ceildiv(tcd_img->x1-tcd_img->x0, tcd_img->comps[compno].dx);
        int i, j;
        for (j=tilec->y0; j<tilec->y1; j++) {
            for (i=tilec->x0; i<tilec->x1; i++) {
                int v;
	#ifdef ENCODE_ONLY
		int v_back;
		v_back =tilec->data[i-tilec->x0+(j-tilec->y0)*tw];
	#endif

                if (tcd_tcp->tccps[compno].qmfbid==1) {
                    v=tilec->data[i-tilec->x0+(j-tilec->y0)*tw];
                } else if (tcd_tcp->tccps[compno].qmfbid==0) {
                    v=tilec->data[i-tilec->x0+(j-tilec->y0)*tw]>>13;
                }
	#ifdef ENCODE_ONLY
			*((double*)matlab_ex_data + i +j*taille_x) = (((double) (v_back))/((double) (1<<13))) ;
	#endif
                v+=adjust;
                tcd_img->comps[compno].data[i+j*w]=int_clamp(v, min, max);

            }
        }

    }
#ifndef TIME_H_NOT_AVAILABLE
    time5=clock()-time5;

    time6=clock()-time6;
#endif

    /*     printf("tile decoding times:\n"); */
    /*     printf("tier 2:    %d.%.3d s\n", time1/CLOCKS_PER_SEC, (time1%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("tier 1:    %d.%.3d s\n", time2/CLOCKS_PER_SEC, (time2%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("dwt:       %d.%.3d s\n", time3/CLOCKS_PER_SEC, (time3%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("mct:       %d.%.3d s\n", time4/CLOCKS_PER_SEC, (time4%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("tile->img: %d.%.3d s\n", time5/CLOCKS_PER_SEC, (time5%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */
    /*     printf("total:     %d.%.3d s\n", time6/CLOCKS_PER_SEC, (time6%CLOCKS_PER_SEC)*1000/CLOCKS_PER_SEC); */

    if (eof) {
        longjmp(j2k_error, 1);
    }
    return l;
}