]> Creatis software - gdcm.git/blob - src/gdcmopenjpeg/libopenjpeg/pi.c
Add IsParsable() method, to avoid confusion with IsReadable()
[gdcm.git] / src / gdcmopenjpeg / libopenjpeg / pi.c
1 /*
2  * Copyright (c) 2001-2002, David Janssens
3  * Copyright (c) 2003-2004, Yannick Verschueren
4  * Copyright (c) 2003-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 "pi.h"
30 #include "int.h"
31 #include <stdlib.h>
32 #include <stdio.h>
33
34 /* <summary>
35  * Create a packet iterator.
36  * </summary> */
37 pi_iterator_t *pi_create(j2k_image_t * img, j2k_cp_t * cp, int tileno)
38 {
39   int p, q, i;
40   int compno, resno, pino;
41   int maxres = 0;
42   pi_iterator_t *pi;
43   j2k_tcp_t *tcp;
44   j2k_tccp_t *tccp;
45
46   tcp = &cp->tcps[tileno];
47   pi = (pi_iterator_t *) malloc((tcp->numpocs + 1) *
48             sizeof(pi_iterator_t));
49
50   for (pino = 0; pino < tcp->numpocs + 1; pino++) {   /* change */
51     p = tileno % cp->tw;
52     q = tileno / cp->tw;
53
54     pi[pino].tx0 = int_max(cp->tx0 + p * cp->tdx, img->x0);
55     pi[pino].ty0 = int_max(cp->ty0 + q * cp->tdy, img->y0);
56     pi[pino].tx1 = int_min(cp->tx0 + (p + 1) * cp->tdx, img->x1);
57     pi[pino].ty1 = int_min(cp->ty0 + (q + 1) * cp->tdy, img->y1);
58     pi[pino].numcomps = img->numcomps;
59     pi[pino].comps =
60       (pi_comp_t *) malloc(img->numcomps * sizeof(pi_comp_t));
61
62     for (compno = 0; compno < pi->numcomps; compno++) {
63       int tcx0, tcy0, tcx1, tcy1;
64       pi_comp_t *comp = &pi[pino].comps[compno];
65       tccp = &tcp->tccps[compno];
66       comp->dx = img->comps[compno].dx;
67       comp->dy = img->comps[compno].dy;
68       comp->numresolutions = tccp->numresolutions;
69       comp->resolutions =
70    (pi_resolution_t *) malloc(comp->numresolutions *
71                sizeof(pi_resolution_t));
72       tcx0 = int_ceildiv(pi->tx0, comp->dx);
73       tcy0 = int_ceildiv(pi->ty0, comp->dy);
74       tcx1 = int_ceildiv(pi->tx1, comp->dx);
75       tcy1 = int_ceildiv(pi->ty1, comp->dy);
76       if (comp->numresolutions > maxres) {
77    maxres = comp->numresolutions;
78       }
79       for (resno = 0; resno < comp->numresolutions; resno++) {
80    int levelno;
81    int rx0, ry0, rx1, ry1;
82    int px0, py0, px1, py1;
83    pi_resolution_t *res = &comp->resolutions[resno];
84    if (tccp->csty & J2K_CCP_CSTY_PRT) {
85      res->pdx = tccp->prcw[resno];
86      res->pdy = tccp->prch[resno];
87    } else {
88      res->pdx = 15;
89      res->pdy = 15;
90    }
91    levelno = comp->numresolutions - 1 - resno;
92    rx0 = int_ceildivpow2(tcx0, levelno);
93    ry0 = int_ceildivpow2(tcy0, levelno);
94    rx1 = int_ceildivpow2(tcx1, levelno);
95    ry1 = int_ceildivpow2(tcy1, levelno);
96    px0 = int_floordivpow2(rx0, res->pdx) << res->pdx;
97    py0 = int_floordivpow2(ry0, res->pdy) << res->pdy;
98    px1 = int_ceildivpow2(rx1, res->pdx) << res->pdx;
99    py1 = int_ceildivpow2(ry1, res->pdy) << res->pdy;
100    res->pw = (rx0==rx1)?0:((px1 - px0) >> res->pdx); /*Mod Antonin : sizebug1*/
101    res->ph = (ry0==ry1)?0:((py1 - py0) >> res->pdy); /*Mod Antonin : sizebug1*/
102       }
103     }
104
105     tccp = &tcp->tccps[0];
106     pi[pino].step_p = 1;
107     pi[pino].step_c = 100 * pi[pino].step_p;
108     pi[pino].step_r = img->numcomps * pi[pino].step_c;
109     pi[pino].step_l = maxres * pi[pino].step_r;
110
111     if (pino == 0) {
112       pi[pino].include =
113    (short int *) malloc(img->numcomps * maxres *
114               tcp->numlayers * 100 * sizeof(short int));
115       for (i = 0; i < img->numcomps * maxres * tcp->numlayers * 100; i++)
116    pi[pino].include[i] = 0;
117     }
118     /* pi[pino].include=(short int*)calloc(img->numcomps*maxres*tcp->numlayers*1000,sizeof(short int)); */
119     else
120       pi[pino].include = pi[pino - 1].include;
121
122     if (tcp->POC == 0) {
123       pi[pino].first = 1;
124       pi[pino].poc.resno0 = 0;
125       pi[pino].poc.compno0 = 0;
126       pi[pino].poc.layno1 = tcp->numlayers;
127       pi[pino].poc.resno1 = maxres;
128       pi[pino].poc.compno1 = img->numcomps;
129       pi[pino].poc.prg = tcp->prg;
130     } else {
131       pi[pino].first = 1;
132       pi[pino].poc.resno0 = tcp->pocs[pino].resno0;
133       pi[pino].poc.compno0 = tcp->pocs[pino].compno0;
134       pi[pino].poc.layno1 = tcp->pocs[pino].layno1;
135       pi[pino].poc.resno1 = tcp->pocs[pino].resno1;
136       pi[pino].poc.compno1 = tcp->pocs[pino].compno1;
137       pi[pino].poc.prg = tcp->pocs[pino].prg;
138     }
139   }
140   return pi;
141 }
142
143 /* <summary>
144  * Get next packet in layer-resolution-component-precinct order.
145  * 
146  * pi: packet iterator to modify
147  * </summary> */
148 int pi_next_lrcp(pi_iterator_t * pi)
149 {
150   pi_comp_t *comp;
151   pi_resolution_t *res;
152
153   if (!pi->first) {
154     comp = &pi->comps[pi->compno];
155     res = &comp->resolutions[pi->resno];
156     goto skip;
157   } else {
158     pi->first = 0;
159   }
160   for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
161     for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1;
162     pi->resno++) {
163       for (pi->compno = pi->poc.compno0;
164       pi->compno < pi->poc.compno1; pi->compno++) {
165    comp = &pi->comps[pi->compno];
166    if (pi->resno >= comp->numresolutions) {
167      continue;
168    }
169    res = &comp->resolutions[pi->resno];
170    for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++) {
171      if (!pi->
172          include[pi->layno * pi->step_l +
173             pi->resno * pi->step_r +
174             pi->compno * pi->step_c + pi->precno * pi->step_p]) {
175        pi->include[pi->layno * pi->step_l +
176          pi->resno * pi->step_r +
177          pi->compno * pi->step_c +
178          pi->precno * pi->step_p] = 1;
179        return 1;
180      }
181    skip:;
182    }
183       }
184     }
185   }
186   return 0;
187 }
188
189 /* <summary>
190  * Get next packet in resolution-layer-component-precinct order.
191  *
192  * pi: packet iterator to modify
193  * </summary> */
194 int pi_next_rlcp(pi_iterator_t * pi)
195 {
196   pi_comp_t *comp;
197   pi_resolution_t *res;
198   if (!pi->first) {
199     comp = &pi->comps[pi->compno];
200     res = &comp->resolutions[pi->resno];
201     goto skip;
202   } else {
203     pi->first = 0;
204   }
205   for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
206     for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
207       for (pi->compno = pi->poc.compno0;
208       pi->compno < pi->poc.compno1; pi->compno++) {
209    comp = &pi->comps[pi->compno];
210    if (pi->resno >= comp->numresolutions) {
211      continue;
212    }
213    res = &comp->resolutions[pi->resno];
214    for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++) {
215      if (!pi->
216          include[pi->layno * pi->step_l +
217             pi->resno * pi->step_r +
218             pi->compno * pi->step_c + pi->precno * pi->step_p]) {
219        pi->include[pi->layno * pi->step_l +
220          pi->resno * pi->step_r +
221          pi->compno * pi->step_c +
222          pi->precno * pi->step_p] = 1;
223        return 1;
224      }
225    skip:;
226    }
227       }
228     }
229   }
230   return 0;
231 }
232
233 /* <summary>
234  * Get next packet in resolution-precinct-component-layer order.
235  *
236  * pi: packet iterator to modify
237  * </summary> */
238 int pi_next_rpcl(pi_iterator_t * pi)
239 {
240   pi_comp_t *comp;
241   pi_resolution_t *res;
242   if (!pi->first) {
243     goto skip;
244   } else {
245     int compno, resno;
246     pi->first = 0;
247     pi->dx = 0;
248     pi->dy = 0;
249     for (compno = 0; compno < pi->numcomps; compno++) {
250       comp = &pi->comps[compno];
251       for (resno = 0; resno < comp->numresolutions; resno++) {
252    int dx, dy;
253    res = &comp->resolutions[resno];
254    dx = comp->dx *
255      (1 << (res->pdx + comp->numresolutions - 1 - resno));
256    dy = comp->dy *
257      (1 << (res->pdy + comp->numresolutions - 1 - resno));
258    pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
259    pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
260       }
261     }
262   }
263   for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
264     for (pi->y = pi->ty0; pi->y < pi->ty1;
265     pi->y += pi->dy - (pi->y % pi->dy)) {
266       for (pi->x = pi->tx0; pi->x < pi->tx1;
267       pi->x += pi->dx - (pi->x % pi->dx)) {
268    for (pi->compno = pi->poc.compno0;
269         pi->compno < pi->poc.compno1; pi->compno++) {
270      int levelno;
271      int trx0, try0;
272      int trx1, try1;/* Add antonin pcrl*/
273      int rpx, rpy;
274      int prci, prcj;
275      comp = &pi->comps[pi->compno];
276      if (pi->resno >= comp->numresolutions) {
277        continue;
278      }
279      res = &comp->resolutions[pi->resno];
280      levelno = comp->numresolutions - 1 - pi->resno;
281      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
282      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
283      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
284      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
285      rpx = res->pdx + levelno;
286      rpy = res->pdy + levelno;
287      if (!
288          (pi->x % (comp->dx << rpx) == 0
289           || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
290        continue;
291      }
292      if (!
293          (pi->y % (comp->dy << rpy) == 0
294           || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
295        continue;
296      }
297
298      /*Add Antonin : sizebug1*/
299      if ((res->pw==0)||(res->pw==0)) continue;
300      /*ddA*/
301
302      /*Add Antonin : pcrl*/
303      if ((trx0==trx1)||(try0==try1)) continue;
304      /*ddA*/
305
306      prci =
307        int_floordivpow2(int_ceildiv
308               (pi->x, comp->dx << levelno),
309               res->pdx) - int_floordivpow2(trx0, res->pdx);
310      prcj =
311        int_floordivpow2(int_ceildiv
312               (pi->y, comp->dy << levelno),
313               res->pdy) - int_floordivpow2(try0, res->pdy);
314      pi->precno = prci + prcj * res->pw;
315      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
316        if (!pi->
317       include[pi->layno * pi->step_l +
318          pi->resno * pi->step_r +
319          pi->compno * pi->step_c +
320          pi->precno * pi->step_p]) {
321          pi->include[pi->layno * pi->step_l +
322            pi->resno * pi->step_r +
323            pi->compno * pi->step_c +
324            pi->precno * pi->step_p] = 1;
325          return 1;
326        }
327      skip:;
328      }
329    }
330       }
331     }
332   }
333   return 0;
334 }
335
336 /* <summary>
337  * Get next packet in precinct-component-resolution-layer order.
338  *
339  * pi: packet iterator to modify
340  * </summary> */
341 int pi_next_pcrl(pi_iterator_t * pi)
342 {
343   pi_comp_t *comp;
344   pi_resolution_t *res;
345   if (!pi->first) {
346     comp = &pi->comps[pi->compno];
347     goto skip;
348   } else {
349     int compno, resno;
350     pi->first = 0;
351     pi->dx = 0;
352     pi->dy = 0;
353     for (compno = 0; compno < pi->numcomps; compno++) {
354       comp = &pi->comps[compno];
355       for (resno = 0; resno < comp->numresolutions; resno++) {
356    int dx, dy;
357    res = &comp->resolutions[resno];
358    dx = comp->dx *
359      (1 << (res->pdx + comp->numresolutions - 1 - resno));
360    dy = comp->dy *
361      (1 << (res->pdy + comp->numresolutions - 1 - resno));
362    pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
363    pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
364       }
365     }
366   }
367   for (pi->y = pi->ty0; pi->y < pi->ty1;
368   pi->y += pi->dy - (pi->y % pi->dy)) {
369     for (pi->x = pi->tx0; pi->x < pi->tx1;
370     pi->x += pi->dx - (pi->x % pi->dx)) {
371       for (pi->compno = pi->poc.compno0;
372       pi->compno < pi->poc.compno1; pi->compno++) {
373    comp = &pi->comps[pi->compno];
374    for (pi->resno = pi->poc.resno0;
375    pi->resno < int_min(pi->poc.resno1,
376      comp->numresolutions); pi->resno++) {
377      int levelno;
378      int trx0, try0;
379      int trx1, try1;/* Add antonin pcrl*/
380      int rpx, rpy;
381      int prci, prcj;
382      res = &comp->resolutions[pi->resno];
383      levelno = comp->numresolutions - 1 - pi->resno;
384      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
385      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
386      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
387      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
388      rpx = res->pdx + levelno;
389      rpy = res->pdy + levelno;
390      if (!
391        (pi->x % (comp->dx << rpx) == 0
392        || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
393        continue;
394      }
395      if (!
396        (pi->y % (comp->dy << rpy) == 0
397        || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
398        continue;
399      }
400      
401      /*Add Antonin : sizebug1*/
402      if ((res->pw==0)||(res->pw==0)) continue;
403      /*ddA*/
404
405      /*Add Antonin : pcrl*/
406      if ((trx0==trx1)||(try0==try1)) continue;
407      /*ddA*/
408      
409      prci =
410        int_floordivpow2(int_ceildiv
411        (pi->x, comp->dx << levelno),
412        res->pdx) - int_floordivpow2(trx0, res->pdx);
413      prcj =
414        int_floordivpow2(int_ceildiv
415        (pi->y, comp->dy << levelno),
416        res->pdy) - int_floordivpow2(try0, res->pdy);
417      pi->precno = prci + prcj * res->pw;
418      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
419        if (!pi->
420          include[pi->layno * pi->step_l +
421          pi->resno * pi->step_r +
422          pi->compno * pi->step_c +
423          pi->precno * pi->step_p]) {
424          pi->include[pi->layno * pi->step_l +
425       pi->resno * pi->step_r +
426       pi->compno * pi->step_c +
427       pi->precno * pi->step_p] = 1;
428          return 1;
429        }
430 skip:;
431      }
432    }
433       }
434     }
435   }
436   return 0;
437 }
438
439 /* <summary>
440  * Get next packet in component-precinct-resolution-layer order.
441  *
442  * pi: packet iterator to modify
443  * </summary> */
444 int pi_next_cprl(pi_iterator_t * pi)
445 {
446   pi_comp_t *comp;
447   pi_resolution_t *res;
448   if (!pi->first) {
449     comp = &pi->comps[pi->compno];
450     goto skip;
451   } else {
452     pi->first = 0;
453   }
454   for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1;
455        pi->compno++) {
456     int resno;
457     comp = &pi->comps[pi->compno];
458     pi->dx = 0;
459     pi->dy = 0;
460     for (resno = 0; resno < comp->numresolutions; resno++) {
461       int dx, dy;
462       res = &comp->resolutions[resno];
463       dx = comp->dx * (1 << (res->pdx + comp->numresolutions - 1 - resno));
464       dy = comp->dy * (1 << (res->pdy + comp->numresolutions - 1 - resno));
465       pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
466       pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
467     }
468     for (pi->y = pi->ty0; pi->y < pi->ty1;
469     pi->y += pi->dy - (pi->y % pi->dy)) {
470       for (pi->x = pi->tx0; pi->x < pi->tx1;
471       pi->x += pi->dx - (pi->x % pi->dx)) {
472    for (pi->resno = pi->poc.resno0;
473         pi->resno < int_min(pi->poc.resno1,
474              comp->numresolutions); pi->resno++) {
475      int levelno;
476      int trx0, try0;
477      int trx1, try1;/* Add antonin pcrl*/
478      int rpx, rpy;
479      int prci, prcj;
480      res = &comp->resolutions[pi->resno];
481      levelno = comp->numresolutions - 1 - pi->resno;
482      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
483      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
484      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
485      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
486      rpx = res->pdx + levelno;
487      rpy = res->pdy + levelno;
488      if (!
489          (pi->x % (comp->dx << rpx) == 0
490           || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
491        continue;
492      }
493      if (!
494          (pi->y % (comp->dy << rpy) == 0
495           || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
496        continue;
497      }
498
499      /*Add Antonin : sizebug1*/
500      if ((res->pw==0)||(res->pw==0)) continue;
501      /*ddA*/
502
503      /*Add Antonin : pcrl*/
504      if ((trx0==trx1)||(try0==try1)) continue;
505      /*ddA*/
506
507      prci =
508        int_floordivpow2(int_ceildiv
509               (pi->x, comp->dx << levelno),
510               res->pdx) - int_floordivpow2(trx0, res->pdx);
511      prcj =
512        int_floordivpow2(int_ceildiv
513               (pi->y, comp->dy << levelno),
514               res->pdy) - int_floordivpow2(try0, res->pdy);
515      pi->precno = prci + prcj * res->pw;
516      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
517        if (!pi->
518       include[pi->layno * pi->step_l +
519          pi->resno * pi->step_r +
520          pi->compno * pi->step_c +
521          pi->precno * pi->step_p]) {
522          pi->include[pi->layno * pi->step_l +
523            pi->resno * pi->step_r +
524            pi->compno * pi->step_c +
525            pi->precno * pi->step_p] = 1;
526          return 1;
527        }
528      skip:;
529      }
530    }
531       }
532     }
533   }
534   return 0;
535 }
536
537 /* <summary>
538  * Get next packet.
539   *
540   * pi: packet iterator to modify
541   * </summary> */
542 int pi_next(pi_iterator_t * pi)
543 {
544   switch (pi->poc.prg) {
545   case 0:
546     return pi_next_lrcp(pi);
547   case 1:
548     return pi_next_rlcp(pi);
549   case 2:
550     return pi_next_rpcl(pi);
551   case 3:
552     return pi_next_pcrl(pi);
553   case 4:
554     return pi_next_cprl(pi);
555   }
556   return 0;
557 }