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