bool tcd_rateallocate(opj_tcd_t *tcd, unsigned char *dest,int len,opj_codestream_info_t *cstr_info)
{
    int compno, resno, bandno, precno, cblkno, passno, layno;
    double min, max;
    double cumdisto[100];    /* fixed_quality */
    const double K = 1;        /* 1.1; fixed_quality */
    double maxSE = 0;

    opj_cp_t *cp = tcd->cp;
    opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
    opj_tcp_t *tcd_tcp = tcd->tcp;

    min = DBL_MAX;
    max = 0;
    
    tcd_tile->numpix = 0;        /* fixed_quality */
   
    for (compno = 0; compno < tcd_tile->numcomps; compno++) {
        opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
        tilec->numpix = 0;

        for (resno = 0; resno < tilec->numresolutions; resno++) {
            opj_tcd_resolution_t *res = &tilec->resolutions[resno];

            for (bandno = 0; bandno < res->numbands; bandno++) {
                opj_tcd_band_t *band = &res->bands[bandno];

                for (precno = 0; precno < res->pw * res->ph; precno++) {
                    opj_tcd_precinct_t *prc = &band->precincts[precno];

                    for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
                        opj_tcd_cblk_enc_t *cblk = &prc->cblks.enc[cblkno];

                        for (passno = 0; passno < cblk->totalpasses; passno++) {
                            opj_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; }
                        } /* passno */
                       
                        /* fixed_quality */
                        tcd_tile->numpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));
                        tilec->numpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));
                    } /* cbklno */
                } /* precno */
            } /* bandno */
        } /* resno */
       
        maxSE += (((double)(1 << tcd->image->comps[compno].prec) - 1.0)
            * ((double)(1 << tcd->image->comps[compno].prec) -1.0))
            * ((double)(tilec->numpix));
    } /* compno */
   
    /* index file */
    if(cstr_info) {
        opj_tile_info_t *tile_info = &cstr_info->tile[tcd->tcd_tileno];
        tile_info->numpix = tcd_tile->numpix;
        tile_info->distotile = tcd_tile->distotile;
        tile_info->thresh = (double *) opj_malloc(tcd_tcp->numlayers * sizeof(double));
    }
   
    for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
        double lo = min;
        double hi = max;
        int success = 0;
        int maxlen = tcd_tcp->rates[layno] ? int_min(((int) ceil(tcd_tcp->rates[layno])), len) : len;
        double goodthresh = 0;
        double stable_thresh = 0;
        int i;
        double distotarget;        /* fixed_quality */
       
        /* fixed_quality */
        distotarget = tcd_tile->distotile - ((K * maxSE) / pow((float)10, tcd_tcp->distoratio[layno] / 10));
       
        /* Don't try to find an optimal threshold but rather take everything not included yet, if
          -r xx,yy,zz,0   (disto_alloc == 1 and rates == 0)
          -q xx,yy,zz,0      (fixed_quality == 1 and distoratio == 0)
          ==> possible to have some lossy layers and the last layer for sure lossless */
        if ( ((cp->disto_alloc==1) && (tcd_tcp->rates[layno]>0)) || ((cp->fixed_quality==1) && (tcd_tcp->distoratio[layno]>0))) {
            opj_t2_t *t2 = t2_create(tcd->cinfo, tcd->image, cp);
            double thresh = 0;

            for (i = 0; i < 32; i++) {
                int l = 0;
                double distoachieved = 0;    /* fixed_quality */
                thresh = (lo + hi) / 2;
               
                tcd_makelayer(tcd, layno, thresh, 0);
               
                if (cp->fixed_quality) {    /* fixed_quality */
                    if(cp->cinema){
                        l = t2_encode_packets(t2,tcd->tcd_tileno, tcd_tile, layno + 1, dest, maxlen, cstr_info,tcd->cur_tp_num,tcd->tp_pos,tcd->cur_pino,THRESH_CALC, tcd->cur_totnum_tp);
                        if (l == -999) {
                            lo = thresh;
                            continue;
                        }else{
                               distoachieved =    layno == 0 ?
                                tcd_tile->distolayer[0]    : cumdisto[layno - 1] + tcd_tile->distolayer[layno];
                            if (distoachieved < distotarget) {
                                hi=thresh;
                                stable_thresh = thresh;
                                continue;
                            }else{
                                lo=thresh;
                            }
                        }
                    }else{
                        distoachieved =    (layno == 0) ?
                            tcd_tile->distolayer[0]    : (cumdisto[layno - 1] + tcd_tile->distolayer[layno]);
                        if (distoachieved < distotarget) {
                            hi = thresh;
                            stable_thresh = thresh;
                            continue;
                        }
                        lo = thresh;
                    }
                } else {
                    l = t2_encode_packets(t2, tcd->tcd_tileno, tcd_tile, layno + 1, dest, maxlen, cstr_info,tcd->cur_tp_num,tcd->tp_pos,tcd->cur_pino,THRESH_CALC, tcd->cur_totnum_tp);
                    /* TODO: what to do with l ??? seek / tell ??? */
                    /* opj_event_msg(tcd->cinfo, EVT_INFO, "rate alloc: len=%d, max=%d\n", l, maxlen); */
                    if (l == -999) {
                        lo = thresh;
                        continue;
                    }
                    hi = thresh;
                    stable_thresh = thresh;
                }
            }
            success = 1;
            goodthresh = stable_thresh == 0? thresh : stable_thresh;
            t2_destroy(t2);
        } else {
            success = 1;
            goodthresh = min;
        }
       
        if (!success) {
            return false;
        }
       
        if(cstr_info) {    /* Threshold for Marcela Index */
            cstr_info->tile[tcd->tcd_tileno].thresh[layno] = goodthresh;
        }
        tcd_makelayer(tcd, layno, goodthresh, 1);
       
        /* fixed_quality */
        cumdisto[layno] = (layno == 0) ? tcd_tile->distolayer[0] : (cumdisto[layno - 1] + tcd_tile->distolayer[layno]);   
    }

    return true;
}