]> Creatis software - gdcm.git/blob - src/gdcmopenjpeg/libopenjpeg/tcd.c
abd0d620d89a7442a3ba14bd61a6aa88a585d29c
[gdcm.git] / src / gdcmopenjpeg / libopenjpeg / tcd.c
1 /*
2  * Copyright (c) 2001-2003, David Janssens
3  * Copyright (c) 2002-2003, Yannick Verschueren
4  * Copyright (c) 2003-2005, Francois Devaux and Antonin Descampe
5  * Copyright (c) 2005, HervĂ© Drolon, FreeImage Team
6  * Copyright (c) 2002-2005, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
7  * All rights reserved.
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions
11  * are met:
12  * 1. Redistributions of source code must retain the above copyright
13  *    notice, this list of conditions and the following disclaimer.
14  * 2. Redistributions in binary form must reproduce the above copyright
15  *    notice, this list of conditions and the following disclaimer in the
16  *    documentation and/or other materials provided with the distribution.
17  *
18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28  * POSSIBILITY OF SUCH DAMAGE.
29  */
30
31 #include "opj_includes.h"
32
33 void tcd_dump(FILE *fd, opj_tcd_t *tcd, opj_tcd_image_t * img) {
34   int tileno, compno, resno, bandno, precno, cblkno;
35
36   fprintf(fd, "image {\n");
37   fprintf(fd, "  tw=%d, th=%d x0=%d x1=%d y0=%d y1=%d\n", 
38     img->tw, img->th, tcd->image->x0, tcd->image->x1, tcd->image->y0, tcd->image->y1);
39
40   for (tileno = 0; tileno < img->th * img->tw; tileno++) {
41     opj_tcd_tile_t *tile = &tcd->tcd_image->tiles[tileno];
42     fprintf(fd, "  tile {\n");
43     fprintf(fd, "    x0=%d, y0=%d, x1=%d, y1=%d, numcomps=%d\n",
44       tile->x0, tile->y0, tile->x1, tile->y1, tile->numcomps);
45     for (compno = 0; compno < tile->numcomps; compno++) {
46       opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
47       fprintf(fd, "    tilec {\n");
48       fprintf(fd,
49         "      x0=%d, y0=%d, x1=%d, y1=%d, numresolutions=%d\n",
50         tilec->x0, tilec->y0, tilec->x1, tilec->y1, tilec->numresolutions);
51       for (resno = 0; resno < tilec->numresolutions; resno++) {
52         opj_tcd_resolution_t *res = &tilec->resolutions[resno];
53         fprintf(fd, "\n   res {\n");
54         fprintf(fd,
55           "          x0=%d, y0=%d, x1=%d, y1=%d, pw=%d, ph=%d, numbands=%d\n",
56           res->x0, res->y0, res->x1, res->y1, res->pw, res->ph, res->numbands);
57         for (bandno = 0; bandno < res->numbands; bandno++) {
58           opj_tcd_band_t *band = &res->bands[bandno];
59           fprintf(fd, "        band {\n");
60           fprintf(fd,
61             "          x0=%d, y0=%d, x1=%d, y1=%d, stepsize=%f, numbps=%d\n",
62             band->x0, band->y0, band->x1, band->y1, band->stepsize, band->numbps);
63           for (precno = 0; precno < res->pw * res->ph; precno++) {
64             opj_tcd_precinct_t *prec = &band->precincts[precno];
65             fprintf(fd, "          prec {\n");
66             fprintf(fd,
67               "            x0=%d, y0=%d, x1=%d, y1=%d, cw=%d, ch=%d\n",
68               prec->x0, prec->y0, prec->x1, prec->y1, prec->cw, prec->ch);
69             for (cblkno = 0; cblkno < prec->cw * prec->ch; cblkno++) {
70               opj_tcd_cblk_t *cblk = &prec->cblks[cblkno];
71               fprintf(fd, "            cblk {\n");
72               fprintf(fd,
73                 "              x0=%d, y0=%d, x1=%d, y1=%d\n",
74                 cblk->x0, cblk->y0, cblk->x1, cblk->y1);
75               fprintf(fd, "            }\n");
76             }
77             fprintf(fd, "          }\n");
78           }
79           fprintf(fd, "        }\n");
80         }
81         fprintf(fd, "      }\n");
82       }
83       fprintf(fd, "    }\n");
84     }
85     fprintf(fd, "  }\n");
86   }
87   fprintf(fd, "}\n");
88 }
89
90 /* ----------------------------------------------------------------------- */
91
92 /**
93 Create a new TCD handle
94 */
95 opj_tcd_t* tcd_create(opj_common_ptr cinfo) {
96   /* create the tcd structure */
97   opj_tcd_t *tcd = (opj_tcd_t*)opj_malloc(sizeof(opj_tcd_t));
98   if(!tcd) return NULL;
99   tcd->cinfo = cinfo;
100   tcd->tcd_image = (opj_tcd_image_t*)opj_malloc(sizeof(opj_tcd_image_t));
101   if(!tcd->tcd_image) {
102     opj_free(tcd);
103     return NULL;
104   }
105
106   return tcd;
107 }
108
109 /**
110 Destroy a previously created TCD handle
111 */
112 void tcd_destroy(opj_tcd_t *tcd) {
113   if(tcd) {
114     opj_free(tcd->tcd_image);
115     opj_free(tcd);
116   }
117 }
118
119 /* ----------------------------------------------------------------------- */
120
121 void tcd_malloc_encode(opj_tcd_t *tcd, opj_image_t * image, opj_cp_t * cp, int curtileno) {
122   int tileno, compno, resno, bandno, precno, cblkno;
123
124   opj_tcd_tile_t *tile = NULL;    /* pointer to tcd->tile */
125   opj_tcd_tilecomp_t *tilec = NULL;  /* pointer to tcd->tilec */
126   opj_tcd_resolution_t *res = NULL;  /* pointer to tcd->res */
127   opj_tcd_band_t *band = NULL;    /* pointer to tcd->band */
128   opj_tcd_precinct_t *prc = NULL;    /* pointer to tcd->prc */
129   opj_tcd_cblk_t *cblk = NULL;    /* pointer to tcd->cblk */
130
131   tcd->image = image;
132   tcd->cp = cp;
133   tcd->tcd_image->tw = cp->tw;
134   tcd->tcd_image->th = cp->th;
135   tcd->tcd_image->tiles = (opj_tcd_tile_t *) opj_malloc(sizeof(opj_tcd_tile_t));
136   
137   for (tileno = 0; tileno < 1; tileno++) {
138     opj_tcp_t *tcp = &cp->tcps[curtileno];
139     int j;
140
141     /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
142     int p = curtileno % cp->tw;  /* si numerotation matricielle .. */
143     int q = curtileno / cp->tw;  /* .. coordonnees de la tile (q,p) q pour ligne et p pour colonne */
144
145     /* opj_tcd_tile_t *tile=&tcd->tcd_image->tiles[tileno]; */
146     tcd->tile = tcd->tcd_image->tiles;
147     tile = tcd->tile;
148
149     /* 4 borders of the tile rescale on the image if necessary */
150     tile->x0 = int_max(cp->tx0 + p * cp->tdx, image->x0);
151     tile->y0 = int_max(cp->ty0 + q * cp->tdy, image->y0);
152     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, image->x1);
153     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, image->y1);
154     tile->numcomps = image->numcomps;
155     /* tile->PPT=image->PPT;  */
156
157     /* Modification of the RATE >> */
158     for (j = 0; j < tcp->numlayers; j++) {
159       tcp->rates[j] = tcp->rates[j] ? 
160         int_ceildiv(tile->numcomps 
161           * (tile->x1 - tile->x0) 
162           * (tile->y1 - tile->y0) 
163           * image->comps[0].prec, 
164           (tcp->rates[j] * 8 * image->comps[0].dx * image->comps[0].dy)) 
165           : 0;
166
167       if (tcp->rates[j]) {
168         if (j && tcp->rates[j] < tcp->rates[j - 1] + 10) {
169           tcp->rates[j] = tcp->rates[j - 1] + 20;
170         } else {
171           if (!j && tcp->rates[j] < 30)
172             tcp->rates[j] = 30;
173         }
174       }
175     }
176     /* << Modification of the RATE */
177     
178     tile->comps = (opj_tcd_tilecomp_t *) opj_malloc(image->numcomps * sizeof(opj_tcd_tilecomp_t));
179     for (compno = 0; compno < tile->numcomps; compno++) {
180       opj_tccp_t *tccp = &tcp->tccps[compno];
181
182       /* opj_tcd_tilecomp_t *tilec=&tile->comps[compno]; */
183       tcd->tilec = &tile->comps[compno];
184       tilec = tcd->tilec;
185
186       /* border of each tile component (global) */
187       tilec->x0 = int_ceildiv(tile->x0, image->comps[compno].dx);
188       tilec->y0 = int_ceildiv(tile->y0, image->comps[compno].dy);
189       tilec->x1 = int_ceildiv(tile->x1, image->comps[compno].dx);
190       tilec->y1 = int_ceildiv(tile->y1, image->comps[compno].dy);
191       
192       tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0) * sizeof(int));
193       tilec->numresolutions = tccp->numresolutions;
194
195       tilec->resolutions = (opj_tcd_resolution_t *) opj_malloc(tilec->numresolutions * sizeof(opj_tcd_resolution_t));
196       
197       for (resno = 0; resno < tilec->numresolutions; resno++) {
198         int pdx, pdy;
199         int levelno = tilec->numresolutions - 1 - resno;
200         int tlprcxstart, tlprcystart, brprcxend, brprcyend;
201         int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
202         int cbgwidthexpn, cbgheightexpn;
203         int cblkwidthexpn, cblkheightexpn;
204
205         /* opj_tcd_resolution_t *res=&tilec->resolutions[resno]; */
206         tcd->res = &tilec->resolutions[resno];
207         res = tcd->res;
208         
209         /* border for each resolution level (global) */
210         res->x0 = int_ceildivpow2(tilec->x0, levelno);
211         res->y0 = int_ceildivpow2(tilec->y0, levelno);
212         res->x1 = int_ceildivpow2(tilec->x1, levelno);
213         res->y1 = int_ceildivpow2(tilec->y1, levelno);
214         
215         res->numbands = resno == 0 ? 1 : 3;
216         /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
217         if (tccp->csty & J2K_CCP_CSTY_PRT) {
218           pdx = tccp->prcw[resno];
219           pdy = tccp->prch[resno];
220         } else {
221           pdx = 15;
222           pdy = 15;
223         }
224         /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
225         tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
226         tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
227         
228         brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
229         brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
230         
231         res->pw = (brprcxend - tlprcxstart) >> pdx;
232         res->ph = (brprcyend - tlprcystart) >> pdy;
233         
234         if (resno == 0) {
235           tlcbgxstart = tlprcxstart;
236           tlcbgystart = tlprcystart;
237           brcbgxend = brprcxend;
238           brcbgyend = brprcyend;
239           cbgwidthexpn = pdx;
240           cbgheightexpn = pdy;
241         } else {
242           tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
243           tlcbgystart = int_ceildivpow2(tlprcystart, 1);
244           brcbgxend = int_ceildivpow2(brprcxend, 1);
245           brcbgyend = int_ceildivpow2(brprcyend, 1);
246           cbgwidthexpn = pdx - 1;
247           cbgheightexpn = pdy - 1;
248         }
249         
250         cblkwidthexpn = int_min(tccp->cblkw, cbgwidthexpn);
251         cblkheightexpn = int_min(tccp->cblkh, cbgheightexpn);
252         
253         for (bandno = 0; bandno < res->numbands; bandno++) {
254           int x0b, y0b, i;
255           int gain, numbps;
256           opj_stepsize_t *ss = NULL;
257
258           tcd->band = &res->bands[bandno];
259           band = tcd->band;
260
261           band->bandno = resno == 0 ? 0 : bandno + 1;
262           x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
263           y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
264           
265           if (band->bandno == 0) {
266             /* band border (global) */
267             band->x0 = int_ceildivpow2(tilec->x0, levelno);
268             band->y0 = int_ceildivpow2(tilec->y0, levelno);
269             band->x1 = int_ceildivpow2(tilec->x1, levelno);
270             band->y1 = int_ceildivpow2(tilec->y1, levelno);
271           } else {
272             /* band border (global) */
273             band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelno) * x0b, levelno + 1);
274             band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelno) * y0b, levelno + 1);
275             band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelno) * x0b, levelno + 1);
276             band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelno) * y0b, levelno + 1);
277           }
278           
279           ss = &tccp->stepsizes[resno == 0 ? 0 : 3 * (resno - 1) + bandno + 1];
280           gain = tccp->qmfbid == 0 ? dwt_getgain_real(band->bandno) : dwt_getgain(band->bandno);          
281           numbps = image->comps[compno].prec + gain;
282           
283           band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0, numbps - ss->expn));
284           band->numbps = ss->expn + tccp->numgbits - 1;  /* WHY -1 ? */
285           
286           band->precincts = (opj_tcd_precinct_t *) opj_malloc(3 * res->pw * res->ph * sizeof(opj_tcd_precinct_t));
287           
288           for (i = 0; i < res->pw * res->ph * 3; i++) {
289             band->precincts[i].imsbtree = NULL;
290             band->precincts[i].incltree = NULL;
291           }
292           
293           for (precno = 0; precno < res->pw * res->ph; precno++) {
294             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
295
296             int cbgxstart = tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
297             int cbgystart = tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
298             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
299             int cbgyend = cbgystart + (1 << cbgheightexpn);
300
301             /* opj_tcd_precinct_t *prc=&band->precincts[precno]; */
302             tcd->prc = &band->precincts[precno];
303             prc = tcd->prc;
304
305             /* precinct size (global) */
306             prc->x0 = int_max(cbgxstart, band->x0);
307             prc->y0 = int_max(cbgystart, band->y0);
308             prc->x1 = int_min(cbgxend, band->x1);
309             prc->y1 = int_min(cbgyend, band->y1);
310
311             tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
312             tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
313             brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
314             brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
315             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
316             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
317
318             prc->cblks = (opj_tcd_cblk_t *) opj_malloc((prc->cw * prc->ch) * sizeof(opj_tcd_cblk_t));
319             prc->incltree = tgt_create(prc->cw, prc->ch);
320             prc->imsbtree = tgt_create(prc->cw, prc->ch);
321             
322             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
323               int cblkxstart = tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
324               int cblkystart = tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
325               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
326               int cblkyend = cblkystart + (1 << cblkheightexpn);
327               
328               tcd->cblk = &prc->cblks[cblkno];
329               cblk = tcd->cblk;
330
331               /* code-block size (global) */
332               cblk->x0 = int_max(cblkxstart, prc->x0);
333               cblk->y0 = int_max(cblkystart, prc->y0);
334               cblk->x1 = int_min(cblkxend, prc->x1);
335               cblk->y1 = int_min(cblkyend, prc->y1);
336             }
337           }
338         }
339       }
340     }
341   }
342   
343   /* tcd_dump(stdout, tcd, &tcd->tcd_image); */
344 }
345
346 void tcd_free_encode(opj_tcd_t *tcd) {
347   int tileno, compno, resno, bandno, precno;
348
349   opj_tcd_tile_t *tile = NULL;    /* pointer to tcd->tile    */
350   opj_tcd_tilecomp_t *tilec = NULL;  /* pointer to tcd->tilec  */
351   opj_tcd_resolution_t *res = NULL;  /* pointer to tcd->res    */
352   opj_tcd_band_t *band = NULL;    /* pointer to tcd->band    */
353   opj_tcd_precinct_t *prc = NULL;    /* pointer to tcd->prc    */
354
355   for (tileno = 0; tileno < 1; tileno++) {
356     tcd->tile = tcd->tcd_image->tiles;
357     tile = tcd->tile;
358
359     for (compno = 0; compno < tile->numcomps; compno++) {
360       tcd->tilec = &tile->comps[compno];
361       tilec = tcd->tilec;
362
363       for (resno = 0; resno < tilec->numresolutions; resno++) {
364         tcd->res = &tilec->resolutions[resno];
365         res = tcd->res;
366
367         for (bandno = 0; bandno < res->numbands; bandno++) {
368           tcd->band = &res->bands[bandno];
369           band = tcd->band;
370
371           for (precno = 0; precno < res->pw * res->ph; precno++) {
372             tcd->prc = &band->precincts[precno];
373             prc = tcd->prc;
374
375             if (prc->incltree != NULL) {
376               tgt_destroy(prc->incltree);
377               prc->incltree = NULL;
378             }
379             if (prc->imsbtree != NULL) {
380               tgt_destroy(prc->imsbtree);  
381               prc->imsbtree = NULL;
382             }
383             opj_free(prc->cblks);
384             prc->cblks = NULL;
385           } /* for (precno */
386           opj_free(band->precincts);
387           band->precincts = NULL;
388         } /* for (bandno */
389       } /* for (resno */
390       opj_free(tilec->resolutions);
391       tilec->resolutions = NULL;
392     } /* for (compno */
393     opj_free(tile->comps);
394     tile->comps = NULL;
395   } /* for (tileno */
396   opj_free(tcd->tcd_image->tiles);
397   tcd->tcd_image->tiles = NULL;
398 }
399
400 void tcd_init_encode(opj_tcd_t *tcd, opj_image_t * image, opj_cp_t * cp, int curtileno) {
401   int tileno, compno, resno, bandno, precno, cblkno;
402
403   opj_tcd_tile_t *tile = NULL;    /* pointer to tcd->tile */
404   opj_tcd_tilecomp_t *tilec = NULL;  /* pointer to tcd->tilec */
405   opj_tcd_resolution_t *res = NULL;  /* pointer to tcd->res */
406   opj_tcd_band_t *band = NULL;    /* pointer to tcd->band */
407   opj_tcd_precinct_t *prc = NULL;    /* pointer to tcd->prc */
408   opj_tcd_cblk_t *cblk = NULL;    /* pointer to tcd->cblk */
409
410   for (tileno = 0; tileno < 1; tileno++) {
411     opj_tcp_t *tcp = &cp->tcps[curtileno];
412     int j;
413     /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
414     int p = curtileno % cp->tw;
415     int q = curtileno / cp->tw;
416
417     tcd->tile = tcd->tcd_image->tiles;
418     tile = tcd->tile;
419     
420     /* 4 borders of the tile rescale on the image if necessary */
421     tile->x0 = int_max(cp->tx0 + p * cp->tdx, image->x0);
422     tile->y0 = int_max(cp->ty0 + q * cp->tdy, image->y0);
423     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, image->x1);
424     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, image->y1);
425     
426     tile->numcomps = image->numcomps;
427     /* tile->PPT=image->PPT; */
428
429     /* Modification of the RATE >> */
430     for (j = 0; j < tcp->numlayers; j++) {
431       tcp->rates[j] = tcp->rates[j] ? 
432         int_ceildiv(tile->numcomps 
433         * (tile->x1 - tile->x0) 
434         * (tile->y1 - tile->y0) 
435         * image->comps[0].prec, 
436         (tcp->rates[j] * 8 * image->comps[0].dx * image->comps[0].dy)) 
437         : 0;
438
439       if (tcp->rates[j]) {
440         if (j && tcp->rates[j] < tcp->rates[j - 1] + 10) {
441           tcp->rates[j] = tcp->rates[j - 1] + 20;
442         } else {
443           if (!j && tcp->rates[j] < 30)
444             tcp->rates[j] = 30;
445         }
446       }
447     }
448     /* << Modification of the RATE */
449
450     /* tile->comps=(opj_tcd_tilecomp_t*)opj_realloc(tile->comps,image->numcomps*sizeof(opj_tcd_tilecomp_t)); */
451     for (compno = 0; compno < tile->numcomps; compno++) {
452       opj_tccp_t *tccp = &tcp->tccps[compno];
453       
454       tcd->tilec = &tile->comps[compno];
455       tilec = tcd->tilec;
456
457       /* border of each tile component (global) */
458       tilec->x0 = int_ceildiv(tile->x0, image->comps[compno].dx);
459       tilec->y0 = int_ceildiv(tile->y0, image->comps[compno].dy);
460       tilec->x1 = int_ceildiv(tile->x1, image->comps[compno].dx);
461       tilec->y1 = int_ceildiv(tile->y1, image->comps[compno].dy);
462       
463       tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0) * sizeof(int));
464       tilec->numresolutions = tccp->numresolutions;
465       /* tilec->resolutions=(opj_tcd_resolution_t*)opj_realloc(tilec->resolutions,tilec->numresolutions*sizeof(opj_tcd_resolution_t)); */
466       for (resno = 0; resno < tilec->numresolutions; resno++) {
467         int pdx, pdy;
468
469         int levelno = tilec->numresolutions - 1 - resno;
470         int tlprcxstart, tlprcystart, brprcxend, brprcyend;
471         int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
472         int cbgwidthexpn, cbgheightexpn;
473         int cblkwidthexpn, cblkheightexpn;
474         
475         tcd->res = &tilec->resolutions[resno];
476         res = tcd->res;
477
478         /* border for each resolution level (global) */
479         res->x0 = int_ceildivpow2(tilec->x0, levelno);
480         res->y0 = int_ceildivpow2(tilec->y0, levelno);
481         res->x1 = int_ceildivpow2(tilec->x1, levelno);
482         res->y1 = int_ceildivpow2(tilec->y1, levelno);  
483         res->numbands = resno == 0 ? 1 : 3;
484
485         /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
486         if (tccp->csty & J2K_CCP_CSTY_PRT) {
487           pdx = tccp->prcw[resno];
488           pdy = tccp->prch[resno];
489         } else {
490           pdx = 15;
491           pdy = 15;
492         }
493         /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
494         tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
495         tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
496         brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
497         brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
498         
499         res->pw = (brprcxend - tlprcxstart) >> pdx;
500         res->ph = (brprcyend - tlprcystart) >> pdy;
501         
502         if (resno == 0) {
503           tlcbgxstart = tlprcxstart;
504           tlcbgystart = tlprcystart;
505           brcbgxend = brprcxend;
506           brcbgyend = brprcyend;
507           cbgwidthexpn = pdx;
508           cbgheightexpn = pdy;
509         } else {
510           tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
511           tlcbgystart = int_ceildivpow2(tlprcystart, 1);
512           brcbgxend = int_ceildivpow2(brprcxend, 1);
513           brcbgyend = int_ceildivpow2(brprcyend, 1);
514           cbgwidthexpn = pdx - 1;
515           cbgheightexpn = pdy - 1;
516         }
517         
518         cblkwidthexpn = int_min(tccp->cblkw, cbgwidthexpn);
519         cblkheightexpn = int_min(tccp->cblkh, cbgheightexpn);
520         
521         for (bandno = 0; bandno < res->numbands; bandno++) {
522           int x0b, y0b;
523           int gain, numbps;
524           opj_stepsize_t *ss = NULL;
525
526           tcd->band = &res->bands[bandno];
527           band = tcd->band;
528
529           band->bandno = resno == 0 ? 0 : bandno + 1;
530           x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
531           y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
532           
533           if (band->bandno == 0) {
534             /* band border */
535             band->x0 = int_ceildivpow2(tilec->x0, levelno);
536             band->y0 = int_ceildivpow2(tilec->y0, levelno);
537             band->x1 = int_ceildivpow2(tilec->x1, levelno);
538             band->y1 = int_ceildivpow2(tilec->y1, levelno);
539           } else {
540             band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelno) * x0b, levelno + 1);
541             band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelno) * y0b, levelno + 1);
542             band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelno) * x0b, levelno + 1);
543             band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelno) * y0b, levelno + 1);
544           }
545           
546           ss = &tccp->stepsizes[resno == 0 ? 0 : 3 * (resno - 1) + bandno + 1];
547           gain = tccp->qmfbid == 0 ? dwt_getgain_real(band->bandno) : dwt_getgain(band->bandno);
548           numbps = image->comps[compno].prec + gain;
549           band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0, numbps - ss->expn));
550           band->numbps = ss->expn + tccp->numgbits - 1;  /* WHY -1 ? */
551           
552           for (precno = 0; precno < res->pw * res->ph; precno++) {
553             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
554
555             int cbgxstart = tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
556             int cbgystart = tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
557             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
558             int cbgyend = cbgystart + (1 << cbgheightexpn);
559             
560             tcd->prc = &band->precincts[precno];
561             prc = tcd->prc;
562
563             /* precinct size (global) */
564             prc->x0 = int_max(cbgxstart, band->x0);
565             prc->y0 = int_max(cbgystart, band->y0);
566             prc->x1 = int_min(cbgxend, band->x1);
567             prc->y1 = int_min(cbgyend, band->y1);
568
569             tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
570             tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
571             brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
572             brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
573             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
574             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
575
576             opj_free(prc->cblks);
577             prc->cblks = (opj_tcd_cblk_t *) opj_malloc(prc->cw * prc->ch * sizeof(opj_tcd_cblk_t));
578
579             if (prc->incltree != NULL) {
580               tgt_destroy(prc->incltree);
581             }
582             if (prc->imsbtree != NULL) {
583               tgt_destroy(prc->imsbtree);
584             }
585             
586             prc->incltree = tgt_create(prc->cw, prc->ch);
587             prc->imsbtree = tgt_create(prc->cw, prc->ch);
588
589             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
590               int cblkxstart = tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
591               int cblkystart = tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
592               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
593               int cblkyend = cblkystart + (1 << cblkheightexpn);
594
595               tcd->cblk = &prc->cblks[cblkno];
596               cblk = tcd->cblk;
597               
598               /* code-block size (global) */
599               cblk->x0 = int_max(cblkxstart, prc->x0);
600               cblk->y0 = int_max(cblkystart, prc->y0);
601               cblk->x1 = int_min(cblkxend, prc->x1);
602               cblk->y1 = int_min(cblkyend, prc->y1);
603             }
604           } /* precno */
605         } /* bandno */
606       } /* resno */
607     } /* compno */
608   } /* tileno */
609
610   /* tcd_dump(stdout, tcd, &tcd->tcd_image); */
611 }
612
613 void tcd_malloc_decode(opj_tcd_t *tcd, opj_image_t * image, opj_cp_t * cp) {
614   int tileno, compno, resno, bandno, precno, cblkno, i, j, p, q;
615   unsigned int x0 = 0, y0 = 0, x1 = 0, y1 = 0, w, h;
616
617   tcd->image = image;
618   tcd->cp = cp;
619   tcd->tcd_image->tw = cp->tw;
620   tcd->tcd_image->th = cp->th;
621   tcd->tcd_image->tiles = (opj_tcd_tile_t *) opj_malloc(cp->tw * cp->th * sizeof(opj_tcd_tile_t));
622   
623   for (i = 0; i < cp->tileno_size; i++) {
624     opj_tcp_t *tcp = &(cp->tcps[cp->tileno[i]]);
625     opj_tcd_tile_t *tile = &(tcd->tcd_image->tiles[cp->tileno[i]]);
626   
627     /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
628     tileno = cp->tileno[i];
629     p = tileno % cp->tw;  /* si numerotation matricielle .. */
630     q = tileno / cp->tw;  /* .. coordonnees de la tile (q,p) q pour ligne et p pour colonne */
631
632     /* 4 borders of the tile rescale on the image if necessary */
633     tile->x0 = int_max(cp->tx0 + p * cp->tdx, image->x0);
634     tile->y0 = int_max(cp->ty0 + q * cp->tdy, image->y0);
635     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, image->x1);
636     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, image->y1);
637     
638     tile->numcomps = image->numcomps;
639     tile->comps = (opj_tcd_tilecomp_t *) opj_malloc(image->numcomps * sizeof(opj_tcd_tilecomp_t));
640     for (compno = 0; compno < tile->numcomps; compno++) {
641       opj_tccp_t *tccp = &tcp->tccps[compno];
642       opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
643
644       /* border of each tile component (global) */
645       tilec->x0 = int_ceildiv(tile->x0, image->comps[compno].dx);
646       tilec->y0 = int_ceildiv(tile->y0, image->comps[compno].dy);
647       tilec->x1 = int_ceildiv(tile->x1, image->comps[compno].dx);
648       tilec->y1 = int_ceildiv(tile->y1, image->comps[compno].dy);
649       
650       tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0) * sizeof(int));
651       tilec->numresolutions = tccp->numresolutions;
652       tilec->resolutions = (opj_tcd_resolution_t *) opj_malloc(tilec->numresolutions * sizeof(opj_tcd_resolution_t));
653
654       for (resno = 0; resno < tilec->numresolutions; resno++) {
655         int pdx, pdy;
656         int levelno = tilec->numresolutions - 1 - resno;
657         int tlprcxstart, tlprcystart, brprcxend, brprcyend;
658         int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
659         int cbgwidthexpn, cbgheightexpn;
660         int cblkwidthexpn, cblkheightexpn;
661
662         opj_tcd_resolution_t *res = &tilec->resolutions[resno];
663         
664         /* border for each resolution level (global) */
665         res->x0 = int_ceildivpow2(tilec->x0, levelno);
666         res->y0 = int_ceildivpow2(tilec->y0, levelno);
667         res->x1 = int_ceildivpow2(tilec->x1, levelno);
668         res->y1 = int_ceildivpow2(tilec->y1, levelno);
669         res->numbands = resno == 0 ? 1 : 3;
670         
671         /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
672         if (tccp->csty & J2K_CCP_CSTY_PRT) {
673           pdx = tccp->prcw[resno];
674           pdy = tccp->prch[resno];
675         } else {
676           pdx = 15;
677           pdy = 15;
678         }
679         
680         /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
681         tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
682         tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
683         brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
684         brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
685
686         res->pw = (res->x0 == res->x1) ? 0 : ((brprcxend - tlprcxstart) >> pdx);
687         res->ph = (res->y0 == res->y1) ? 0 : ((brprcyend - tlprcystart) >> pdy);
688         
689         if (resno == 0) {
690           tlcbgxstart = tlprcxstart;
691           tlcbgystart = tlprcystart;
692           brcbgxend = brprcxend;
693           brcbgyend = brprcyend;
694           cbgwidthexpn = pdx;
695           cbgheightexpn = pdy;
696         } else {
697           tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
698           tlcbgystart = int_ceildivpow2(tlprcystart, 1);
699           brcbgxend = int_ceildivpow2(brprcxend, 1);
700           brcbgyend = int_ceildivpow2(brprcyend, 1);
701           cbgwidthexpn = pdx - 1;
702           cbgheightexpn = pdy - 1;
703         }
704         
705         cblkwidthexpn = int_min(tccp->cblkw, cbgwidthexpn);
706         cblkheightexpn = int_min(tccp->cblkh, cbgheightexpn);
707         
708         for (bandno = 0; bandno < res->numbands; bandno++) {
709           int x0b, y0b;
710           int gain, numbps;
711           opj_stepsize_t *ss = NULL;
712
713           opj_tcd_band_t *band = &res->bands[bandno];
714           band->bandno = resno == 0 ? 0 : bandno + 1;
715           x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
716           y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
717           
718           if (band->bandno == 0) {
719             /* band border (global) */
720             band->x0 = int_ceildivpow2(tilec->x0, levelno);
721             band->y0 = int_ceildivpow2(tilec->y0, levelno);
722             band->x1 = int_ceildivpow2(tilec->x1, levelno);
723             band->y1 = int_ceildivpow2(tilec->y1, levelno);
724           } else {
725             /* band border (global) */
726             band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelno) * x0b, levelno + 1);
727             band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelno) * y0b, levelno + 1);
728             band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelno) * x0b, levelno + 1);
729             band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelno) * y0b, levelno + 1);
730           }
731           
732           ss = &tccp->stepsizes[resno == 0 ? 0 : 3 * (resno - 1) + bandno + 1];
733           gain = tccp->qmfbid == 0 ? dwt_getgain_real(band->bandno) : dwt_getgain(band->bandno);
734           numbps = image->comps[compno].prec + gain;
735           band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0, numbps - ss->expn));
736           band->numbps = ss->expn + tccp->numgbits - 1;  /* WHY -1 ? */
737           
738           band->precincts = (opj_tcd_precinct_t *) opj_malloc(res->pw * res->ph * sizeof(opj_tcd_precinct_t));
739           
740           for (precno = 0; precno < res->pw * res->ph; precno++) {
741             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
742             int cbgxstart = tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
743             int cbgystart = tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
744             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
745             int cbgyend = cbgystart + (1 << cbgheightexpn);
746
747             opj_tcd_precinct_t *prc = &band->precincts[precno];
748             /* precinct size (global) */
749             prc->x0 = int_max(cbgxstart, band->x0);
750             prc->y0 = int_max(cbgystart, band->y0);
751             prc->x1 = int_min(cbgxend, band->x1);
752             prc->y1 = int_min(cbgyend, band->y1);
753             
754             tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
755             tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
756             brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
757             brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
758             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
759             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
760             
761             prc->cblks = (opj_tcd_cblk_t *) opj_malloc(prc->cw * prc->ch * sizeof(opj_tcd_cblk_t));
762             
763             prc->incltree = tgt_create(prc->cw, prc->ch);
764             prc->imsbtree = tgt_create(prc->cw, prc->ch);
765             
766             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
767               int cblkxstart = tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
768               int cblkystart = tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
769               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
770               int cblkyend = cblkystart + (1 << cblkheightexpn);          
771               
772               /* code-block size (global) */
773               opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
774               cblk->x0 = int_max(cblkxstart, prc->x0);
775               cblk->y0 = int_max(cblkystart, prc->y0);
776               cblk->x1 = int_min(cblkxend, prc->x1);
777               cblk->y1 = int_min(cblkyend, prc->y1);
778             }
779           } /* precno */
780         } /* bandno */
781       } /* resno */
782     } /* compno */
783   } /* i = 0..cp->tileno_size */
784
785   /* tcd_dump(stdout, tcd, &tcd->tcd_image); */
786
787   /* 
788   Allocate place to store the decoded data = final image
789   Place limited by the tile really present in the codestream 
790   */
791   
792   for (i = 0; i < image->numcomps; i++) {
793     for (j = 0; j < cp->tileno_size; j++) {
794       tileno = cp->tileno[j];
795       x0 = j == 0 ? tcd->tcd_image->tiles[tileno].comps[i].x0 : int_min(x0,
796         (unsigned int) tcd->tcd_image->tiles[tileno].comps[i].x0);
797       y0 = j == 0 ? tcd->tcd_image->tiles[tileno].comps[i].y0 : int_min(y0,
798         (unsigned int) tcd->tcd_image->tiles[tileno].comps[i].y0);
799       x1 = j == 0 ? tcd->tcd_image->tiles[tileno].comps[i].x1 : int_max(x1,
800         (unsigned int) tcd->tcd_image->tiles[tileno].comps[i].x1);
801       y1 = j == 0 ? tcd->tcd_image->tiles[tileno].comps[i].y1 : int_max(y1, 
802         (unsigned int) tcd->tcd_image->tiles[tileno].comps[i].y1);
803     }
804     
805     w = x1 - x0;
806     h = y1 - y0;
807     
808     image->comps[i].data = (int *) opj_malloc(w * h * sizeof(int));
809     image->comps[i].w = w;
810     image->comps[i].h = h;
811     image->comps[i].x0 = x0;
812     image->comps[i].y0 = y0;
813   }
814 }
815
816 void tcd_makelayer_fixed(opj_tcd_t *tcd, int layno, int final) {
817   int compno, resno, bandno, precno, cblkno;
818   int value;      /*, matrice[tcd_tcp->numlayers][tcd_tile->comps[0].numresolutions][3]; */
819   int matrice[10][10][3];
820   int i, j, k;
821
822   opj_cp_t *cp = tcd->cp;
823   opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
824   opj_tcp_t *tcd_tcp = tcd->tcp;
825
826   /*matrice=(int*)opj_malloc(tcd_tcp->numlayers*tcd_tile->comps[0].numresolutions*3*sizeof(int)); */
827   
828   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
829     opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
830     for (i = 0; i < tcd_tcp->numlayers; i++) {
831       for (j = 0; j < tilec->numresolutions; j++) {
832         for (k = 0; k < 3; k++) {
833           matrice[i][j][k] =
834             (int) (cp->matrice[i * tilec->numresolutions * 3 + j * 3 + k] 
835             * (float) (tcd->image->comps[compno].prec / 16.0));
836         }
837       }
838     }
839         
840     for (resno = 0; resno < tilec->numresolutions; resno++) {
841       opj_tcd_resolution_t *res = &tilec->resolutions[resno];
842       for (bandno = 0; bandno < res->numbands; bandno++) {
843         opj_tcd_band_t *band = &res->bands[bandno];
844         for (precno = 0; precno < res->pw * res->ph; precno++) {
845           opj_tcd_precinct_t *prc = &band->precincts[precno];
846           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
847             opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
848             opj_tcd_layer_t *layer = &cblk->layers[layno];
849             int n;
850             int imsb = tcd->image->comps[compno].prec - cblk->numbps;  /* number of bit-plan equal to zero */
851             /* Correction of the matrix of coefficient to include the IMSB information */
852             if (layno == 0) {
853               value = matrice[layno][resno][bandno];
854               if (imsb >= value) {
855                 value = 0;
856               } else {
857                 value -= imsb;
858               }
859             } else {
860               value =  matrice[layno][resno][bandno] -  matrice[layno - 1][resno][bandno];
861               if (imsb >= matrice[layno - 1][resno][bandno]) {
862                 value -= (imsb - matrice[layno - 1][resno][bandno]);
863                 if (value < 0) {
864                   value = 0;
865                 }
866               }
867             }
868             
869             if (layno == 0) {
870               cblk->numpassesinlayers = 0;
871             }
872             
873             n = cblk->numpassesinlayers;
874             if (cblk->numpassesinlayers == 0) {
875               if (value != 0) {
876                 n = 3 * value - 2 + cblk->numpassesinlayers;
877               } else {
878                 n = cblk->numpassesinlayers;
879               }
880             } else {
881               n = 3 * value + cblk->numpassesinlayers;
882             }
883             
884             layer->numpasses = n - cblk->numpassesinlayers;
885             
886             if (!layer->numpasses)
887               continue;
888             
889             if (cblk->numpassesinlayers == 0) {
890               layer->len = cblk->passes[n - 1].rate;
891               layer->data = cblk->data;
892             } else {
893               layer->len = cblk->passes[n - 1].rate - cblk->passes[cblk->numpassesinlayers - 1].rate;
894               layer->data = cblk->data + cblk->passes[cblk->numpassesinlayers - 1].rate;
895             }
896             if (final)
897               cblk->numpassesinlayers = n;
898           }
899         }
900       }
901     }
902   }
903 }
904
905 void tcd_rateallocate_fixed(opj_tcd_t *tcd) {
906   int layno;
907   for (layno = 0; layno < tcd->tcp->numlayers; layno++) {
908     tcd_makelayer_fixed(tcd, layno, 1);
909   }
910 }
911
912 void tcd_makelayer(opj_tcd_t *tcd, int layno, double thresh, int final) {
913   int compno, resno, bandno, precno, cblkno, passno;
914   
915   opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
916
917   tcd_tile->distolayer[layno] = 0;  /* fixed_quality */
918   
919   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
920     opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
921     for (resno = 0; resno < tilec->numresolutions; resno++) {
922       opj_tcd_resolution_t *res = &tilec->resolutions[resno];
923       for (bandno = 0; bandno < res->numbands; bandno++) {
924         opj_tcd_band_t *band = &res->bands[bandno];
925         for (precno = 0; precno < res->pw * res->ph; precno++) {
926           opj_tcd_precinct_t *prc = &band->precincts[precno];
927           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
928             opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
929             opj_tcd_layer_t *layer = &cblk->layers[layno];
930             
931             int n;
932             if (layno == 0) {
933               cblk->numpassesinlayers = 0;
934             }
935             n = cblk->numpassesinlayers;
936             for (passno = cblk->numpassesinlayers; passno < cblk->totalpasses; passno++) {
937               int dr;
938               double dd;
939               opj_tcd_pass_t *pass = &cblk->passes[passno];
940               if (n == 0) {
941                 dr = pass->rate;
942                 dd = pass->distortiondec;
943               } else {
944                 dr = pass->rate - cblk->passes[n - 1].rate;
945                 dd = pass->distortiondec - cblk->passes[n - 1].distortiondec;
946               }
947               if (!dr) {
948                 if (dd != 0)
949                   n = passno + 1;
950                 continue;
951               }
952               if (dd / dr >= thresh)
953                 n = passno + 1;
954             }
955             layer->numpasses = n - cblk->numpassesinlayers;
956             
957             if (!layer->numpasses) {
958               layer->disto = 0;
959               continue;
960             }
961             if (cblk->numpassesinlayers == 0) {
962               layer->len = cblk->passes[n - 1].rate;
963               layer->data = cblk->data;
964               layer->disto = cblk->passes[n - 1].distortiondec;
965             } else {
966               layer->len = cblk->passes[n - 1].rate -  cblk->passes[cblk->numpassesinlayers - 1].rate;
967               layer->data = cblk->data + cblk->passes[cblk->numpassesinlayers - 1].rate;
968               layer->disto = cblk->passes[n - 1].distortiondec - cblk->passes[cblk->numpassesinlayers - 1].distortiondec;
969             }
970             
971             tcd_tile->distolayer[layno] += layer->disto;  /* fixed_quality */
972             
973             if (final)
974               cblk->numpassesinlayers = n;
975           }
976         }
977       }
978     }
979   }
980 }
981
982 bool tcd_rateallocate(opj_tcd_t *tcd, unsigned char *dest, int len, opj_image_info_t * image_info) {
983   int compno, resno, bandno, precno, cblkno, passno, layno;
984   double min, max;
985   double cumdisto[100];  /* fixed_quality */
986   const double K = 1;    /* 1.1; // fixed_quality */
987   double maxSE = 0;
988
989   opj_cp_t *cp = tcd->cp;
990   opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
991   opj_tcp_t *tcd_tcp = tcd->tcp;
992
993   min = DBL_MAX;
994   max = 0;
995   
996   tcd_tile->nbpix = 0;    /* fixed_quality */
997   
998   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
999     opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1000     tilec->nbpix = 0;
1001     for (resno = 0; resno < tilec->numresolutions; resno++) {
1002       opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1003       for (bandno = 0; bandno < res->numbands; bandno++) {
1004         opj_tcd_band_t *band = &res->bands[bandno];
1005         for (precno = 0; precno < res->pw * res->ph; precno++) {
1006           opj_tcd_precinct_t *prc = &band->precincts[precno];
1007           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
1008             opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1009             for (passno = 0; passno < cblk->totalpasses; passno++) {
1010               opj_tcd_pass_t *pass = &cblk->passes[passno];
1011               int dr;
1012               double dd, rdslope;
1013               if (passno == 0) {
1014                 dr = pass->rate;
1015                 dd = pass->distortiondec;
1016               } else {
1017                 dr = pass->rate - cblk->passes[passno - 1].rate;
1018                 dd = pass->distortiondec - cblk->passes[passno - 1].distortiondec;
1019               }
1020               if (dr == 0) {
1021                 continue;
1022               }
1023               rdslope = dd / dr;
1024               if (rdslope < min) {
1025                 min = rdslope;
1026               }
1027               if (rdslope > max) {
1028                 max = rdslope;
1029               }
1030             } /* passno */
1031             
1032             /* fixed_quality */
1033             tcd_tile->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));
1034             tilec->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));
1035           } /* cbklno */
1036         } /* precno */
1037       } /* bandno */
1038     } /* resno */
1039     
1040     maxSE += (((double)(1 << tcd->image->comps[compno].prec) - 1.0) 
1041       * ((double)(1 << tcd->image->comps[compno].prec) -1.0)) 
1042       * ((double)(tilec->nbpix));
1043   } /* compno */
1044   
1045   /* add antonin index */
1046   if(image_info && image_info->index_on) {
1047     opj_tile_info_t *info_TL = &image_info->tile[tcd->tcd_tileno];
1048     info_TL->nbpix = tcd_tile->nbpix;
1049     info_TL->distotile = tcd_tile->distotile;
1050     info_TL->thresh = (double *) opj_malloc(tcd_tcp->numlayers * sizeof(double));
1051   }
1052   /* dda */
1053   
1054   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
1055     double lo = min;
1056     double hi = max;
1057     int success = 0;
1058     int maxlen = tcd_tcp->rates[layno] ? int_min(tcd_tcp->rates[layno], len) : len;
1059     double goodthresh = 0;
1060     int i;
1061     double distotarget;    /* fixed_quality */
1062     
1063     /* fixed_quality */
1064     distotarget = tcd_tile->distotile - ((K * maxSE) / pow((float)10, tcd_tcp->distoratio[layno] / 10));
1065         
1066     if ((tcd_tcp->rates[layno]) || (cp->disto_alloc==0)) {
1067       opj_t2_t *t2 = t2_create(tcd->cinfo, tcd->image, cp);
1068
1069       for (i = 0; i < 32; i++) {
1070         double thresh = (lo + hi) / 2;
1071         int l = 0;
1072         double distoachieved = 0;  /* fixed_quality */
1073         
1074         tcd_makelayer(tcd, layno, thresh, 0);
1075         
1076         if (cp->fixed_quality) {  /* fixed_quality */
1077           distoachieved =  layno == 0 ? 
1078             tcd_tile->distolayer[0]  : cumdisto[layno - 1] + tcd_tile->distolayer[layno];
1079           if (distoachieved < distotarget) {
1080             hi = thresh;
1081             continue;
1082           }
1083           lo = thresh;
1084         } else {
1085           l = t2_encode_packets(t2, tcd->tcd_tileno, tcd_tile, layno + 1, dest, maxlen, image_info);
1086           /* opj_event_msg(tcd->cinfo, EVT_INFO, "rate alloc: len=%d, max=%d\n", l, maxlen); */
1087           if (l == -999) {
1088             lo = thresh;
1089             continue;
1090           }
1091           hi = thresh;
1092         }
1093         
1094         success = 1;
1095         goodthresh = thresh;
1096       }
1097       t2_destroy(t2);
1098     } else {
1099       success = 1;
1100       goodthresh = min;
1101     }
1102     
1103     if (!success) {
1104       return false;
1105     }
1106     
1107     if(image_info && image_info->index_on) {  /* Threshold for Marcela Index */
1108       image_info->tile[tcd->tcd_tileno].thresh[layno] = goodthresh;
1109     }
1110     tcd_makelayer(tcd, layno, goodthresh, 1);
1111         
1112     /* fixed_quality */
1113     cumdisto[layno] = layno == 0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] + tcd_tile->distolayer[layno];  
1114   }
1115
1116   return true;
1117 }
1118
1119 int tcd_encode_tile(opj_tcd_t *tcd, int tileno, unsigned char *dest, int len, opj_image_info_t * image_info) {
1120   int compno;
1121   int l, i, npck = 0;
1122   double encoding_time;
1123   opj_tcd_tile_t *tile = NULL;
1124   opj_tcp_t *tcd_tcp = NULL;
1125   opj_cp_t *cp = NULL;
1126
1127   opj_tcp_t *tcp = &tcd->cp->tcps[0];
1128   opj_tccp_t *tccp = &tcp->tccps[0];
1129   opj_image_t *image = tcd->image;
1130   
1131   opj_t1_t *t1 = NULL;    /* T1 component */
1132   opj_t2_t *t2 = NULL;    /* T2 component */
1133
1134   tcd->tcd_tileno = tileno;
1135   tcd->tcd_tile = tcd->tcd_image->tiles;
1136   tcd->tcp = &tcd->cp->tcps[tileno];
1137
1138   tile = tcd->tcd_tile;
1139   tcd_tcp = tcd->tcp;
1140   cp = tcd->cp;
1141
1142   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1143   if(image_info && image_info->index_on) {
1144     opj_tcd_tilecomp_t *tilec_idx = &tile->comps[0];  /* based on component 0 */
1145     for (i = 0; i < tilec_idx->numresolutions; i++) {
1146       opj_tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1147
1148       image_info->tile[tileno].pw[i] = res_idx->pw;
1149       image_info->tile[tileno].ph[i] = res_idx->ph;
1150
1151       npck += res_idx->pw * res_idx->ph;
1152
1153       image_info->tile[tileno].pdx[i] = tccp->prcw[i];
1154       image_info->tile[tileno].pdy[i] = tccp->prch[i];
1155     }
1156     image_info->tile[tileno].packet = (opj_packet_info_t *) opj_malloc(image_info->comp * image_info->layer * npck * sizeof(opj_packet_info_t));
1157   }
1158   /* << INDEX */
1159   
1160   /*---------------TILE-------------------*/
1161   encoding_time = opj_clock();  /* time needed to encode a tile */
1162   
1163   for (compno = 0; compno < tile->numcomps; compno++) {
1164     int x, y;
1165
1166     int adjust = image->comps[compno].sgnd ? 0 : 1 << (image->comps[compno].prec - 1);
1167     int offset_x = int_ceildiv(image->x0, image->comps[compno].dx);
1168     int offset_y = int_ceildiv(image->y0, image->comps[compno].dy);
1169     
1170     opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1171     int tw = tilec->x1 - tilec->x0;
1172     int w = int_ceildiv(image->x1 - image->x0, image->comps[compno].dx);
1173
1174     /* extract tile data */
1175
1176     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1177       for (y = tilec->y0; y < tilec->y1; y++) {
1178         /* start of the src tile scanline */
1179         int *data = &image->comps[compno].data[(tilec->x0 - offset_x) + (y - offset_y) * w];
1180         /* start of the dst tile scanline */
1181         int *tile_data = &tilec->data[(y - tilec->y0) * tw];
1182         for (x = tilec->x0; x < tilec->x1; x++) {
1183           *tile_data++ = *data++ - adjust;
1184         }
1185       }
1186     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1187       for (y = tilec->y0; y < tilec->y1; y++) {
1188         /* start of the src tile scanline */
1189         int *data = &image->comps[compno].data[(tilec->x0 - offset_x) + (y - offset_y) * w];
1190         /* start of the dst tile scanline */
1191         int *tile_data = &tilec->data[(y - tilec->y0) * tw];
1192         for (x = tilec->x0; x < tilec->x1; x++) {
1193           *tile_data++ = (*data++ - adjust) << 13;
1194         }
1195       }
1196     }
1197   }
1198   
1199   /*----------------MCT-------------------*/
1200   if (tcd_tcp->mct) {
1201     int samples = (tile->comps[0].x1 - tile->comps[0].x0) * (tile->comps[0].y1 - tile->comps[0].y0);
1202     if (tcd_tcp->tccps[0].qmfbid == 0) {
1203       mct_encode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, samples);
1204     } else {
1205       mct_encode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, samples);
1206     }
1207   }
1208   
1209   /*----------------DWT---------------------*/
1210
1211   for (compno = 0; compno < tile->numcomps; compno++) {
1212     opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1213     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1214       dwt_encode(tilec);
1215     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1216       dwt_encode_real(tilec);
1217     }
1218   }
1219
1220   /*------------------TIER1-----------------*/
1221   t1 = t1_create(tcd->cinfo);
1222   t1_encode_cblks(t1, tile, tcd_tcp);
1223   t1_destroy(t1);
1224   
1225   /*-----------RATE-ALLOCATE------------------*/
1226
1227   /* INDEX */
1228   if(image_info) {
1229     image_info->index_write = 0;
1230   }
1231   if (cp->disto_alloc || cp->fixed_quality) {  /* fixed_quality */
1232     /* Normal Rate/distortion allocation */
1233     tcd_rateallocate(tcd, dest, len, image_info);
1234   } else {
1235     /* Fixed layer allocation */
1236     tcd_rateallocate_fixed(tcd);
1237   }
1238   
1239   /*--------------TIER2------------------*/
1240
1241   /* INDEX */
1242   if(image_info) {
1243     image_info->index_write = 1;
1244   }
1245
1246   t2 = t2_create(tcd->cinfo, image, cp);
1247   l = t2_encode_packets(t2, tileno, tile, tcd_tcp->numlayers, dest, len, image_info);
1248   t2_destroy(t2);
1249   
1250   /*---------------CLEAN-------------------*/
1251
1252   encoding_time = opj_clock() - encoding_time;
1253   opj_event_msg(tcd->cinfo, EVT_INFO, "- tile encoded in %f s\n", encoding_time);
1254   
1255   /* cleaning memory */
1256   for (compno = 0; compno < tile->numcomps; compno++) {
1257     tcd->tilec = &tile->comps[compno];
1258     opj_free(tcd->tilec->data);
1259   }
1260   
1261   return l;
1262 }
1263
1264 bool tcd_decode_tile(opj_tcd_t *tcd, unsigned char *src, int len, int tileno) {
1265   int l;
1266   int compno;
1267   int eof = 0;
1268   double tile_time, t1_time, dwt_time;
1269   opj_tcd_tile_t *tile = NULL;
1270
1271   opj_t1_t *t1 = NULL;    /* T1 component */
1272   opj_t2_t *t2 = NULL;    /* T2 component */
1273   
1274   tcd->tcd_tileno = tileno;
1275   tcd->tcd_tile = &(tcd->tcd_image->tiles[tileno]);
1276   tcd->tcp = &(tcd->cp->tcps[tileno]);
1277   tile = tcd->tcd_tile;
1278   
1279   tile_time = opj_clock();  /* time needed to decode a tile */
1280   opj_event_msg(tcd->cinfo, EVT_INFO, "tile %d of %d\n", tileno + 1, tcd->cp->tw * tcd->cp->th);
1281   
1282   /*--------------TIER2------------------*/
1283   
1284   t2 = t2_create(tcd->cinfo, tcd->image, tcd->cp);
1285   l = t2_decode_packets(t2, src, len, tileno, tile);
1286   t2_destroy(t2);
1287
1288   if (l == -999) {
1289     eof = 1;
1290     opj_event_msg(tcd->cinfo, EVT_ERROR, "tcd_decode: incomplete bistream\n");
1291   }
1292   
1293   /*------------------TIER1-----------------*/
1294   
1295   t1_time = opj_clock();  /* time needed to decode a tile */
1296   t1 = t1_create(tcd->cinfo);
1297   t1_decode_cblks(t1, tile, tcd->tcp);
1298   t1_destroy(t1);
1299   t1_time = opj_clock() - t1_time;
1300   opj_event_msg(tcd->cinfo, EVT_INFO, "- tiers-1 took %f s\n", t1_time);
1301   
1302   /*----------------DWT---------------------*/
1303
1304   dwt_time = opj_clock();  /* time needed to decode a tile */
1305   for (compno = 0; compno < tile->numcomps; compno++) {
1306     opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1307     if (tcd->cp->reduce != 0) {
1308       tcd->image->comps[compno].resno_decoded =
1309         tile->comps[compno].numresolutions - tcd->cp->reduce - 1;
1310     }
1311         
1312     if (tcd->tcp->tccps[compno].qmfbid == 1) {
1313       dwt_decode(tilec, tilec->numresolutions - 1 - tcd->image->comps[compno].resno_decoded);
1314     } else {
1315       dwt_decode_real(tilec, tilec->numresolutions - 1 - tcd->image->comps[compno].resno_decoded);
1316     }
1317
1318     if (tile->comps[compno].numresolutions > 0) {
1319       tcd->image->comps[compno].factor = tile->comps[compno].numresolutions - (tcd->image->comps[compno].resno_decoded + 1);
1320     }
1321   }
1322   dwt_time = opj_clock() - dwt_time;
1323   opj_event_msg(tcd->cinfo, EVT_INFO, "- dwt took %f s\n", dwt_time);
1324   
1325   /*----------------MCT-------------------*/
1326   
1327   if (tcd->tcp->mct) {
1328     if (tcd->tcp->tccps[0].qmfbid == 1) {
1329       mct_decode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, 
1330         (tile->comps[0].x1 - tile->comps[0].x0) * (tile->comps[0].y1 - tile->comps[0].y0));
1331     } else {
1332       mct_decode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data, 
1333         (tile->comps[0].x1 - tile->comps[0].x0) * (tile->comps[0].y1 - tile->comps[0].y0));
1334     }
1335   }
1336   
1337   /*---------------TILE-------------------*/
1338   
1339   for (compno = 0; compno < tile->numcomps; compno++) {
1340     opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1341     opj_tcd_resolution_t *res =  &tilec->resolutions[tcd->image->comps[compno].resno_decoded];
1342     int adjust = tcd->image->comps[compno].sgnd ? 0 : 1 << (tcd->image->comps[compno].prec - 1);
1343     int min = tcd->image->comps[compno].sgnd ? 
1344       -(1 << (tcd->image->comps[compno].prec - 1)) : 0;
1345     int max = tcd->image->comps[compno].sgnd ? 
1346       (1 << (tcd->image->comps[compno].prec - 1)) - 1 : (1 << tcd->image->comps[compno].prec) - 1;
1347     
1348     int tw = tilec->x1 - tilec->x0;
1349     int w = tcd->image->comps[compno].w;
1350     
1351     int i, j;
1352     int offset_x = int_ceildivpow2(tcd->image->comps[compno].x0, tcd->image->comps[compno].factor);
1353     int offset_y = int_ceildivpow2(tcd->image->comps[compno].y0, tcd->image->comps[compno].factor);
1354     
1355     for (j = res->y0; j < res->y1; j++) {
1356       for (i = res->x0; i < res->x1; i++) {
1357         int v;
1358         float tmp = (float)((tilec->data[i - res->x0 + (j - res->y0) * tw]) / 8192.0);
1359
1360         if (tcd->tcp->tccps[compno].qmfbid == 1) {
1361           v = tilec->data[i - res->x0 + (j - res->y0) * tw];
1362         } else {
1363           int tmp2 = ((int) (floor(fabs(tmp)))) + ((int) floor(fabs(tmp*2))%2);
1364           v = ((tmp < 0) ? -tmp2:tmp2);
1365         }
1366         v += adjust;
1367         
1368         tcd->image->comps[compno].data[(i - offset_x) + (j - offset_y) * w] = int_clamp(v, min, max);
1369       }
1370     }
1371   }
1372   
1373   tile_time = opj_clock() - tile_time;  /* time needed to decode a tile */
1374   opj_event_msg(tcd->cinfo, EVT_INFO, "- tile decoded in %f s\n", tile_time);
1375     
1376   for (compno = 0; compno < tile->numcomps; compno++) {
1377     opj_free(tcd->tcd_image->tiles[tileno].comps[compno].data);
1378     tcd->tcd_image->tiles[tileno].comps[compno].data = NULL;
1379   }
1380   
1381   if (eof) {
1382     return false;
1383   }
1384   
1385   return true;
1386 }
1387
1388 void tcd_free_decode(opj_tcd_t *tcd) {
1389   int tileno,compno,resno,bandno,precno;
1390
1391   opj_tcd_image_t *tcd_image = tcd->tcd_image;
1392   
1393   for (tileno = 0; tileno < tcd_image->tw * tcd_image->th; tileno++) {
1394     opj_tcd_tile_t *tile = &tcd_image->tiles[tileno];
1395     for (compno = 0; compno < tile->numcomps; compno++) {
1396       opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1397       for (resno = 0; resno < tilec->numresolutions; resno++) {
1398         opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1399         for (bandno = 0; bandno < res->numbands; bandno++) {
1400           opj_tcd_band_t *band = &res->bands[bandno];
1401           for (precno = 0; precno < res->ph * res->pw; precno++) {
1402             opj_tcd_precinct_t *prec = &band->precincts[precno];
1403             if (prec->cblks != NULL) opj_free(prec->cblks);
1404             if (prec->imsbtree != NULL) tgt_destroy(prec->imsbtree);
1405             if (prec->incltree != NULL) tgt_destroy(prec->incltree);
1406           }
1407           if (band->precincts != NULL) opj_free(band->precincts);
1408         }
1409       }
1410       if (tilec->resolutions != NULL) opj_free(tilec->resolutions);
1411     }
1412     if (tile->comps != NULL) opj_free(tile->comps);
1413   }
1414
1415   if (tcd_image->tiles != NULL) opj_free(tcd_image->tiles);
1416 }
1417