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