src/gdcmjasper/src/libjasper/jpc/jpc_qmfb.c

   1 /*
   2  * Copyright (c) 1999-2000 Image Power, Inc. and the University of
   3  *   British Columbia.
   4  * Copyright (c) 2001-2003 Michael David Adams.
   5  * All rights reserved.
   6  */
   7
   8 /* __START_OF_JASPER_LICENSE__
   9  *
  10  * JasPer License Version 2.0
  11  *
  12  * Copyright (c) 1999-2000 Image Power, Inc.
  13  * Copyright (c) 1999-2000 The University of British Columbia
  14  * Copyright (c) 2001-2003 Michael David Adams
  15  *
  16  * All rights reserved.
  17  *
  18  * Permission is hereby granted, free of charge, to any person (the
  19  * "User") obtaining a copy of this software and associated documentation
  20  * files (the "Software"), to deal in the Software without restriction,
  21  * including without limitation the rights to use, copy, modify, merge,
  22  * publish, distribute, and/or sell copies of the Software, and to permit
  23  * persons to whom the Software is furnished to do so, subject to the
  24  * following conditions:
  25  *
  26  * 1.  The above copyright notices and this permission notice (which
  27  * includes the disclaimer below) shall be included in all copies or
  28  * substantial portions of the Software.
  29  *
  30  * 2.  The name of a copyright holder shall not be used to endorse or
  31  * promote products derived from the Software without specific prior
  32  * written permission.
  33  *
  34  * THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL PART OF THIS
  35  * LICENSE.  NO USE OF THE SOFTWARE IS AUTHORIZED HEREUNDER EXCEPT UNDER
  36  * THIS DISCLAIMER.  THE SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS
  37  * "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
  38  * BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
  39  * PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.  IN NO
  40  * EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL
  41  * INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING
  42  * FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
  43  * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
  44  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.  NO ASSURANCES ARE
  45  * PROVIDED BY THE COPYRIGHT HOLDERS THAT THE SOFTWARE DOES NOT INFRINGE
  46  * THE PATENT OR OTHER INTELLECTUAL PROPERTY RIGHTS OF ANY OTHER ENTITY.
  47  * EACH COPYRIGHT HOLDER DISCLAIMS ANY LIABILITY TO THE USER FOR CLAIMS
  48  * BROUGHT BY ANY OTHER ENTITY BASED ON INFRINGEMENT OF INTELLECTUAL
  49  * PROPERTY RIGHTS OR OTHERWISE.  AS A CONDITION TO EXERCISING THE RIGHTS
  50  * GRANTED HEREUNDER, EACH USER HEREBY ASSUMES SOLE RESPONSIBILITY TO SECURE
  51  * ANY OTHER INTELLECTUAL PROPERTY RIGHTS NEEDED, IF ANY.  THE SOFTWARE
  52  * IS NOT FAULT-TOLERANT AND IS NOT INTENDED FOR USE IN MISSION-CRITICAL
  53  * SYSTEMS, SUCH AS THOSE USED IN THE OPERATION OF NUCLEAR FACILITIES,
  54  * AIRCRAFT NAVIGATION OR COMMUNICATION SYSTEMS, AIR TRAFFIC CONTROL
  55  * SYSTEMS, DIRECT LIFE SUPPORT MACHINES, OR WEAPONS SYSTEMS, IN WHICH
  56  * THE FAILURE OF THE SOFTWARE OR SYSTEM COULD LEAD DIRECTLY TO DEATH,
  57  * PERSONAL INJURY, OR SEVERE PHYSICAL OR ENVIRONMENTAL DAMAGE ("HIGH
  58  * RISK ACTIVITIES").  THE COPYRIGHT HOLDERS SPECIFICALLY DISCLAIM ANY
  59  * EXPRESS OR IMPLIED WARRANTY OF FITNESS FOR HIGH RISK ACTIVITIES.
  60  *
  61  * __END_OF_JASPER_LICENSE__
  62  */
  63
  64 /*
  65  * Quadrature Mirror-Image Filter Bank (QMFB) Library
  66  *
  67  * $Id: jpc_qmfb.c,v 1.1 2005/05/22 18:33:05 malaterre Exp $
  68  */
  69
  70 /******************************************************************************\
  71 * Includes.
  72 \******************************************************************************/
  73
  74 #include <assert.h>
  75
  76 #include "jasper/jas_fix.h"
  77 #include "jasper/jas_malloc.h"
  78 #include "jasper/jas_math.h"
  79
  80 #include "jpc_qmfb.h"
  81 #include "jpc_tsfb.h"
  82 #include "jpc_math.h"
  83
  84 /******************************************************************************\
  85 *
  86 \******************************************************************************/
  87
  88 static jpc_qmfb1d_t *jpc_qmfb1d_create(void);
  89
  90 static int jpc_ft_getnumchans(jpc_qmfb1d_t *qmfb);
  91 static int jpc_ft_getanalfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters);
  92 static int jpc_ft_getsynfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters);
  93 static void jpc_ft_analyze(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x);
  94 static void jpc_ft_synthesize(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x);
  95
  96 static int jpc_ns_getnumchans(jpc_qmfb1d_t *qmfb);
  97 static int jpc_ns_getanalfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters);
  98 static int jpc_ns_getsynfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters);
  99 static void jpc_ns_analyze(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x);
 100 static void jpc_ns_synthesize(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x);
 101
 102 /******************************************************************************\
 103 *
 104 \******************************************************************************/
 105
 106 jpc_qmfb1dops_t jpc_ft_ops = {
 107   jpc_ft_getnumchans,
 108   jpc_ft_getanalfilters,
 109   jpc_ft_getsynfilters,
 110   jpc_ft_analyze,
 111   jpc_ft_synthesize
 112 };
 113
 114 jpc_qmfb1dops_t jpc_ns_ops = {
 115   jpc_ns_getnumchans,
 116   jpc_ns_getanalfilters,
 117   jpc_ns_getsynfilters,
 118   jpc_ns_analyze,
 119   jpc_ns_synthesize
 120 };
 121
 122 /******************************************************************************\
 123 *
 124 \******************************************************************************/
 125
 126 static void jpc_qmfb1d_setup(jpc_fix_t *startptr, int startind, int endind,
 127   int intrastep, jpc_fix_t **lstartptr, int *lstartind, int *lendind,
 128   jpc_fix_t **hstartptr, int *hstartind, int *hendind)
 129 {
 130   *lstartind = JPC_CEILDIVPOW2(startind, 1);
 131   *lendind = JPC_CEILDIVPOW2(endind, 1);
 132   *hstartind = JPC_FLOORDIVPOW2(startind, 1);
 133   *hendind = JPC_FLOORDIVPOW2(endind, 1);
 134   *lstartptr = startptr;
 135   *hstartptr = &startptr[(*lendind - *lstartind) * intrastep];
 136 }
 137
 138 static void jpc_qmfb1d_split(jpc_fix_t *startptr, int startind, int endind,
 139   register int step, jpc_fix_t *lstartptr, int lstartind, int lendind,
 140   jpc_fix_t *hstartptr, int hstartind, int hendind)
 141 {
 142   int bufsize = JPC_CEILDIVPOW2(endind - startind, 2);
 143 #if !defined(HAVE_VLA)
 144 #define QMFB_SPLITBUFSIZE 4096
 145   jpc_fix_t splitbuf[QMFB_SPLITBUFSIZE];
 146 #else
 147   jpc_fix_t splitbuf[bufsize];
 148 #endif
 149   jpc_fix_t *buf = splitbuf;
 150   int llen;
 151   int hlen;
 152   int twostep;
 153   jpc_fix_t *tmpptr;
 154   register jpc_fix_t *ptr;
 155   register jpc_fix_t *hptr;
 156   register jpc_fix_t *lptr;
 157   register int n;
 158   int state;
 159
 160   twostep = step << 1;
 161   llen = lendind - lstartind;
 162   hlen = hendind - hstartind;
 163
 164 #if !defined(HAVE_VLA)
 165   /* Get a buffer. */
 166   if (bufsize > QMFB_SPLITBUFSIZE) {
 167     if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
 168       /* We have no choice but to commit suicide in this case. */
 169       abort();
 170     }
 171   }
 172 #endif
 173
 174   if (hstartind < lstartind) {
 175     /* The first sample in the input signal is to appear
 176       in the highpass subband signal. */
 177     /* Copy the appropriate samples into the lowpass subband
 178       signal, saving any samples destined for the highpass subband
 179       signal as they are overwritten. */
 180     tmpptr = buf;
 181     ptr = &startptr[step];
 182     lptr = lstartptr;
 183     n = llen;
 184     state = 1;
 185     while (n-- > 0) {
 186       if (state) {
 187         *tmpptr = *lptr;
 188         ++tmpptr;
 189       }
 190       *lptr = *ptr;
 191       ptr += twostep;
 192       lptr += step;
 193       state ^= 1;
 194     }
 195     /* Copy the appropriate samples into the highpass subband
 196       signal. */
 197     /* Handle the nonoverwritten samples. */
 198     hptr = &hstartptr[(hlen - 1) * step];
 199     ptr = &startptr[(((llen + hlen - 1) >> 1) << 1) * step];
 200     n = hlen - (tmpptr - buf);
 201     while (n-- > 0) {
 202       *hptr = *ptr;
 203       hptr -= step;
 204       ptr -= twostep;
 205     }
 206     /* Handle the overwritten samples. */
 207     n = tmpptr - buf;
 208     while (n-- > 0) {
 209       --tmpptr;
 210       *hptr = *tmpptr;
 211       hptr -= step;
 212     }
 213   } else {
 214     /* The first sample in the input signal is to appear
 215       in the lowpass subband signal. */
 216     /* Copy the appropriate samples into the lowpass subband
 217       signal, saving any samples for the highpass subband
 218       signal as they are overwritten. */
 219     state = 0;
 220     ptr = startptr;
 221     lptr = lstartptr;
 222     tmpptr = buf;
 223     n = llen;
 224     while (n-- > 0) {
 225       if (state) {
 226         *tmpptr = *lptr;
 227         ++tmpptr;
 228       }
 229       *lptr = *ptr;
 230       ptr += twostep;
 231       lptr += step;
 232       state ^= 1;
 233     }
 234     /* Copy the appropriate samples into the highpass subband
 235       signal. */
 236     /* Handle the nonoverwritten samples. */
 237     ptr = &startptr[((((llen + hlen) >> 1) << 1) - 1) * step];
 238     hptr = &hstartptr[(hlen - 1) * step];
 239     n = hlen - (tmpptr - buf);
 240     while (n-- > 0) {
 241       *hptr = *ptr;
 242       ptr -= twostep;
 243       hptr -= step;
 244     }
 245     /* Handle the overwritten samples. */
 246     n = tmpptr - buf;
 247     while (n-- > 0) {
 248       --tmpptr;
 249       *hptr = *tmpptr;
 250       hptr -= step;
 251     }
 252   }
 253
 254 #if !defined(HAVE_VLA)
 255   /* If the split buffer was allocated on the heap, free this memory. */
 256   if (buf != splitbuf) {
 257     jas_free(buf);
 258   }
 259 #endif
 260 }
 261
 262 static void jpc_qmfb1d_join(jpc_fix_t *startptr, int startind, int endind,
 263   register int step, jpc_fix_t *lstartptr, int lstartind, int lendind,
 264   jpc_fix_t *hstartptr, int hstartind, int hendind)
 265 {
 266   int bufsize = JPC_CEILDIVPOW2(endind - startind, 2);
 267 #if !defined(HAVE_VLA)
 268 #define  QMFB_JOINBUFSIZE  4096
 269   jpc_fix_t joinbuf[QMFB_JOINBUFSIZE];
 270 #else
 271   jpc_fix_t joinbuf[bufsize];
 272 #endif
 273   jpc_fix_t *buf = joinbuf;
 274   int llen;
 275   int hlen;
 276   int twostep;
 277   jpc_fix_t *tmpptr;
 278   register jpc_fix_t *ptr;
 279   register jpc_fix_t *hptr;
 280   register jpc_fix_t *lptr;
 281   register int n;
 282   int state;
 283
 284 #if !defined(HAVE_VLA)
 285   /* Allocate memory for the join buffer from the heap. */
 286   if (bufsize > QMFB_JOINBUFSIZE) {
 287     if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
 288       /* We have no choice but to commit suicide. */
 289       abort();
 290     }
 291   }
 292 #endif
 293
 294   twostep = step << 1;
 295   llen = lendind - lstartind;
 296   hlen = hendind - hstartind;
 297
 298   if (hstartind < lstartind) {
 299     /* The first sample in the highpass subband signal is to
 300       appear first in the output signal. */
 301     /* Copy the appropriate samples into the first phase of the
 302       output signal. */
 303     tmpptr = buf;
 304     hptr = hstartptr;
 305     ptr = startptr;
 306     n = (llen + 1) >> 1;
 307     while (n-- > 0) {
 308       *tmpptr = *ptr;
 309       *ptr = *hptr;
 310       ++tmpptr;
 311       ptr += twostep;
 312       hptr += step;
 313     }
 314     n = hlen - ((llen + 1) >> 1);
 315     while (n-- > 0) {
 316       *ptr = *hptr;
 317       ptr += twostep;
 318       hptr += step;
 319     }
 320     /* Copy the appropriate samples into the second phase of
 321       the output signal. */
 322     ptr -= (lendind > hendind) ? (step) : (step + twostep);
 323     state = !((llen - 1) & 1);
 324     lptr = &lstartptr[(llen - 1) * step];
 325     n = llen;
 326     while (n-- > 0) {
 327       if (state) {
 328         --tmpptr;
 329         *ptr = *tmpptr;
 330       } else {
 331         *ptr = *lptr;
 332       }
 333       lptr -= step;
 334       ptr -= twostep;
 335       state ^= 1;
 336     }
 337   } else {
 338     /* The first sample in the lowpass subband signal is to
 339       appear first in the output signal. */
 340     /* Copy the appropriate samples into the first phase of the
 341       output signal (corresponding to even indexed samples). */
 342     lptr = &lstartptr[(llen - 1) * step];
 343     ptr = &startptr[((llen - 1) << 1) * step];
 344     n = llen >> 1;
 345     tmpptr = buf;
 346     while (n-- > 0) {
 347       *tmpptr = *ptr;
 348       *ptr = *lptr;
 349       ++tmpptr;
 350       ptr -= twostep;
 351       lptr -= step;
 352     }
 353     n = llen - (llen >> 1);
 354     while (n-- > 0) {
 355       *ptr = *lptr;
 356       ptr -= twostep;
 357       lptr -= step;
 358     }
 359     /* Copy the appropriate samples into the second phase of
 360       the output signal (corresponding to odd indexed
 361       samples). */
 362     ptr = &startptr[step];
 363     hptr = hstartptr;
 364     state = !(llen & 1);
 365     n = hlen;
 366     while (n-- > 0) {
 367       if (state) {
 368         --tmpptr;
 369         *ptr = *tmpptr;
 370       } else {
 371         *ptr = *hptr;
 372       }
 373       hptr += step;
 374       ptr += twostep;
 375       state ^= 1;
 376     }
 377   }
 378
 379 #if !defined(HAVE_VLA)
 380   /* If the join buffer was allocated on the heap, free this memory. */
 381   if (buf != joinbuf) {
 382     jas_free(buf);
 383   }
 384 #endif
 385 }
 386
 387 /******************************************************************************\
 388 * Code for 5/3 transform.
 389 \******************************************************************************/
 390
 391 static int jpc_ft_getnumchans(jpc_qmfb1d_t *qmfb)
 392 {
 393   /* Avoid compiler warnings about unused parameters. */
 394   qmfb = 0;
 395
 396   return 2;
 397 }
 398
 399 static int jpc_ft_getanalfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
 400 {
 401   /* Avoid compiler warnings about unused parameters. */
 402   qmfb = 0;
 403   len = 0;
 404   filters = 0;
 405   abort();
 406   return -1;
 407 }
 408
 409 static int jpc_ft_getsynfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
 410 {
 411   jas_seq_t *lf;
 412   jas_seq_t *hf;
 413
 414   /* Avoid compiler warnings about unused parameters. */
 415   qmfb = 0;
 416
 417   lf = 0;
 418   hf = 0;
 419
 420   if (len > 1 || (!len)) {
 421     if (!(lf = jas_seq_create(-1, 2))) {
 422       goto error;
 423     }
 424     jas_seq_set(lf, -1, jpc_dbltofix(0.5));
 425     jas_seq_set(lf, 0, jpc_dbltofix(1.0));
 426     jas_seq_set(lf, 1, jpc_dbltofix(0.5));
 427     if (!(hf = jas_seq_create(-1, 4))) {
 428       goto error;
 429     }
 430     jas_seq_set(hf, -1, jpc_dbltofix(-0.125));
 431     jas_seq_set(hf, 0, jpc_dbltofix(-0.25));
 432     jas_seq_set(hf, 1, jpc_dbltofix(0.75));
 433     jas_seq_set(hf, 2, jpc_dbltofix(-0.25));
 434     jas_seq_set(hf, 3, jpc_dbltofix(-0.125));
 435   } else if (len == 1) {
 436     if (!(lf = jas_seq_create(0, 1))) {
 437       goto error;
 438     }
 439     jas_seq_set(lf, 0, jpc_dbltofix(1.0));
 440     if (!(hf = jas_seq_create(0, 1))) {
 441       goto error;
 442     }
 443     jas_seq_set(hf, 0, jpc_dbltofix(2.0));
 444   } else {
 445     abort();
 446   }
 447
 448   filters[0] = lf;
 449   filters[1] = hf;
 450
 451   return 0;
 452
 453 error:
 454   if (lf) {
 455     jas_seq_destroy(lf);
 456   }
 457   if (hf) {
 458     jas_seq_destroy(hf);
 459   }
 460   return -1;
 461 }
 462
 463 #define  NFT_LIFT0(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, pluseq) \
 464 { \
 465   register jpc_fix_t *lptr = (lstartptr); \
 466   register jpc_fix_t *hptr = (hstartptr); \
 467   register int n = (hendind) - (hstartind); \
 468   if ((hstartind) < (lstartind)) { \
 469     pluseq(*hptr, *lptr); \
 470     hptr += (step); \
 471     --n; \
 472   } \
 473   if ((hendind) >= (lendind)) { \
 474     --n; \
 475   } \
 476   while (n-- > 0) { \
 477     pluseq(*hptr, jpc_fix_asr(jpc_fix_add(*lptr, lptr[(step)]), 1)); \
 478     hptr += (step); \
 479     lptr += (step); \
 480   } \
 481   if ((hendind) >= (lendind)) { \
 482     pluseq(*hptr, *lptr); \
 483   } \
 484 }
 485
 486 #define  NFT_LIFT1(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, pluseq) \
 487 { \
 488   register jpc_fix_t *lptr = (lstartptr); \
 489   register jpc_fix_t *hptr = (hstartptr); \
 490   register int n = (lendind) - (lstartind); \
 491   if ((hstartind) >= (lstartind)) { \
 492     pluseq(*lptr, *hptr); \
 493     lptr += (step); \
 494     --n; \
 495   } \
 496   if ((lendind) > (hendind)) { \
 497     --n; \
 498   } \
 499   while (n-- > 0) { \
 500     pluseq(*lptr, jpc_fix_asr(jpc_fix_add(*hptr, hptr[(step)]), 2)); \
 501     lptr += (step); \
 502     hptr += (step); \
 503   } \
 504   if ((lendind) > (hendind)) { \
 505     pluseq(*lptr, *hptr); \
 506   } \
 507 }
 508
 509 #define  RFT_LIFT0(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, pmeqop) \
 510 { \
 511   register jpc_fix_t *lptr = (lstartptr); \
 512   register jpc_fix_t *hptr = (hstartptr); \
 513   register int n = (hendind) - (hstartind); \
 514   if ((hstartind) < (lstartind)) { \
 515     *hptr pmeqop *lptr; \
 516     hptr += (step); \
 517     --n; \
 518   } \
 519   if ((hendind) >= (lendind)) { \
 520     --n; \
 521   } \
 522   while (n-- > 0) { \
 523     *hptr pmeqop (*lptr + lptr[(step)]) >> 1; \
 524     hptr += (step); \
 525     lptr += (step); \
 526   } \
 527   if ((hendind) >= (lendind)) { \
 528     *hptr pmeqop *lptr; \
 529   } \
 530 }
 531
 532 #define  RFT_LIFT1(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, pmeqop) \
 533 { \
 534   register jpc_fix_t *lptr = (lstartptr); \
 535   register jpc_fix_t *hptr = (hstartptr); \
 536   register int n = (lendind) - (lstartind); \
 537   if ((hstartind) >= (lstartind)) { \
 538     *lptr pmeqop ((*hptr << 1) + 2) >> 2; \
 539     lptr += (step); \
 540     --n; \
 541   } \
 542   if ((lendind) > (hendind)) { \
 543     --n; \
 544   } \
 545   while (n-- > 0) { \
 546     *lptr pmeqop ((*hptr + hptr[(step)]) + 2) >> 2; \
 547     lptr += (step); \
 548     hptr += (step); \
 549   } \
 550   if ((lendind) > (hendind)) { \
 551     *lptr pmeqop ((*hptr << 1) + 2) >> 2; \
 552   } \
 553 }
 554
 555 static void jpc_ft_analyze(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
 556 {
 557   jpc_fix_t *startptr;
 558   int startind;
 559   int endind;
 560   jpc_fix_t *  lstartptr;
 561   int   lstartind;
 562   int   lendind;
 563   jpc_fix_t *  hstartptr;
 564   int   hstartind;
 565   int   hendind;
 566   int interstep;
 567   int intrastep;
 568   int numseq;
 569
 570   /* Avoid compiler warnings about unused parameters. */
 571   qmfb = 0;
 572
 573   if (flags & JPC_QMFB1D_VERT) {
 574     interstep = 1;
 575     intrastep = jas_seq2d_rowstep(x);
 576     numseq = jas_seq2d_width(x);
 577     startind = jas_seq2d_ystart(x);
 578     endind = jas_seq2d_yend(x);
 579   } else {
 580     interstep = jas_seq2d_rowstep(x);
 581     intrastep = 1;
 582     numseq = jas_seq2d_height(x);
 583     startind = jas_seq2d_xstart(x);
 584     endind = jas_seq2d_xend(x);
 585   }
 586
 587   assert(startind < endind);
 588
 589   startptr = jas_seq2d_getref(x, jas_seq2d_xstart(x), jas_seq2d_ystart(x));
 590   if (flags & JPC_QMFB1D_RITIMODE) {
 591     while (numseq-- > 0) {
 592       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 593         &lstartptr, &lstartind, &lendind, &hstartptr,
 594         &hstartind, &hendind);
 595       if (endind - startind > 1) {
 596         jpc_qmfb1d_split(startptr, startind, endind,
 597           intrastep, lstartptr, lstartind, lendind,
 598           hstartptr, hstartind, hendind);
 599         RFT_LIFT0(lstartptr, lstartind, lendind,
 600           hstartptr, hstartind, hendind, intrastep, -=);
 601         RFT_LIFT1(lstartptr, lstartind, lendind,
 602           hstartptr, hstartind, hendind, intrastep, +=);
 603       } else {
 604         if (lstartind == lendind) {
 605           *startptr <<= 1;
 606         }
 607       }
 608       startptr += interstep;
 609     }
 610   } else {
 611     while (numseq-- > 0) {
 612       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 613         &lstartptr, &lstartind, &lendind, &hstartptr,
 614         &hstartind, &hendind);
 615       if (endind - startind > 1) {
 616         jpc_qmfb1d_split(startptr, startind, endind,
 617           intrastep, lstartptr, lstartind, lendind,
 618           hstartptr, hstartind, hendind);
 619         NFT_LIFT0(lstartptr, lstartind, lendind,
 620           hstartptr, hstartind, hendind, intrastep,
 621           jpc_fix_minuseq);
 622         NFT_LIFT1(lstartptr, lstartind, lendind,
 623           hstartptr, hstartind, hendind, intrastep,
 624           jpc_fix_pluseq);
 625       } else {
 626         if (lstartind == lendind) {
 627           *startptr = jpc_fix_asl(*startptr, 1);
 628         }
 629       }
 630       startptr += interstep;
 631     }
 632   }
 633 }
 634
 635 static void jpc_ft_synthesize(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
 636 {
 637   jpc_fix_t *startptr;
 638   int startind;
 639   int endind;
 640   jpc_fix_t *lstartptr;
 641   int lstartind;
 642   int lendind;
 643   jpc_fix_t *hstartptr;
 644   int hstartind;
 645   int hendind;
 646   int interstep;
 647   int intrastep;
 648   int numseq;
 649
 650   /* Avoid compiler warnings about unused parameters. */
 651   qmfb = 0;
 652
 653   if (flags & JPC_QMFB1D_VERT) {
 654     interstep = 1;
 655     intrastep = jas_seq2d_rowstep(x);
 656     numseq = jas_seq2d_width(x);
 657     startind = jas_seq2d_ystart(x);
 658     endind = jas_seq2d_yend(x);
 659   } else {
 660     interstep = jas_seq2d_rowstep(x);
 661     intrastep = 1;
 662     numseq = jas_seq2d_height(x);
 663     startind = jas_seq2d_xstart(x);
 664     endind = jas_seq2d_xend(x);
 665   }
 666
 667   assert(startind < endind);
 668
 669   startptr = jas_seq2d_getref(x, jas_seq2d_xstart(x), jas_seq2d_ystart(x));
 670   if (flags & JPC_QMFB1D_RITIMODE) {
 671     while (numseq-- > 0) {
 672       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 673         &lstartptr, &lstartind, &lendind, &hstartptr,
 674         &hstartind, &hendind);
 675       if (endind - startind > 1) {
 676         RFT_LIFT1(lstartptr, lstartind, lendind,
 677           hstartptr, hstartind, hendind, intrastep, -=);
 678         RFT_LIFT0(lstartptr, lstartind, lendind,
 679           hstartptr, hstartind, hendind, intrastep, +=);
 680         jpc_qmfb1d_join(startptr, startind, endind,
 681           intrastep, lstartptr, lstartind, lendind,
 682           hstartptr, hstartind, hendind);
 683       } else {
 684         if (lstartind == lendind) {
 685           *startptr >>= 1;
 686         }
 687       }
 688       startptr += interstep;
 689     }
 690   } else {
 691     while (numseq-- > 0) {
 692       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 693         &lstartptr, &lstartind, &lendind, &hstartptr,
 694         &hstartind, &hendind);
 695       if (endind - startind > 1) {
 696         NFT_LIFT1(lstartptr, lstartind, lendind,
 697           hstartptr, hstartind, hendind, intrastep,
 698           jpc_fix_minuseq);
 699         NFT_LIFT0(lstartptr, lstartind, lendind,
 700           hstartptr, hstartind, hendind, intrastep,
 701           jpc_fix_pluseq);
 702         jpc_qmfb1d_join(startptr, startind, endind,
 703           intrastep, lstartptr, lstartind, lendind,
 704           hstartptr, hstartind, hendind);
 705       } else {
 706         if (lstartind == lendind) {
 707           *startptr = jpc_fix_asr(*startptr, 1);
 708         }
 709       }
 710       startptr += interstep;
 711     }
 712   }
 713 }
 714
 715 /******************************************************************************\
 716 * Code for 9/7 transform.
 717 \******************************************************************************/
 718
 719 static int jpc_ns_getnumchans(jpc_qmfb1d_t *qmfb)
 720 {
 721   /* Avoid compiler warnings about unused parameters. */
 722   qmfb = 0;
 723
 724   return 2;
 725 }
 726
 727 static int jpc_ns_getanalfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
 728 {
 729   /* Avoid compiler warnings about unused parameters. */
 730   qmfb = 0;
 731   len = 0;
 732   filters = 0;
 733
 734   abort();
 735   return -1;
 736 }
 737
 738 static int jpc_ns_getsynfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
 739 {
 740   jas_seq_t *lf;
 741   jas_seq_t *hf;
 742
 743   /* Avoid compiler warnings about unused parameters. */
 744   qmfb = 0;
 745
 746   lf = 0;
 747   hf = 0;
 748
 749   if (len > 1 || (!len)) {
 750     if (!(lf = jas_seq_create(-3, 4))) {
 751       goto error;
 752     }
 753     jas_seq_set(lf, -3, jpc_dbltofix(-0.09127176311424948));
 754     jas_seq_set(lf, -2, jpc_dbltofix(-0.05754352622849957));
 755     jas_seq_set(lf, -1, jpc_dbltofix(0.5912717631142470));
 756     jas_seq_set(lf, 0, jpc_dbltofix(1.115087052456994));
 757     jas_seq_set(lf, 1, jpc_dbltofix(0.5912717631142470));
 758     jas_seq_set(lf, 2, jpc_dbltofix(-0.05754352622849957));
 759     jas_seq_set(lf, 3, jpc_dbltofix(-0.09127176311424948));
 760     if (!(hf = jas_seq_create(-3, 6))) {
 761       goto error;
 762     }
 763     jas_seq_set(hf, -3, jpc_dbltofix(-0.02674875741080976 * 2.0));
 764     jas_seq_set(hf, -2, jpc_dbltofix(-0.01686411844287495 * 2.0));
 765     jas_seq_set(hf, -1, jpc_dbltofix(0.07822326652898785 * 2.0));
 766     jas_seq_set(hf, 0, jpc_dbltofix(0.2668641184428723 * 2.0));
 767     jas_seq_set(hf, 1, jpc_dbltofix(-0.6029490182363579 * 2.0));
 768     jas_seq_set(hf, 2, jpc_dbltofix(0.2668641184428723 * 2.0));
 769     jas_seq_set(hf, 3, jpc_dbltofix(0.07822326652898785 * 2.0));
 770     jas_seq_set(hf, 4, jpc_dbltofix(-0.01686411844287495 * 2.0));
 771     jas_seq_set(hf, 5, jpc_dbltofix(-0.02674875741080976 * 2.0));
 772   } else if (len == 1) {
 773     if (!(lf = jas_seq_create(0, 1))) {
 774       goto error;
 775     }
 776     jas_seq_set(lf, 0, jpc_dbltofix(1.0));
 777     if (!(hf = jas_seq_create(0, 1))) {
 778       goto error;
 779     }
 780     jas_seq_set(hf, 0, jpc_dbltofix(2.0));
 781   } else {
 782     abort();
 783   }
 784
 785   filters[0] = lf;
 786   filters[1] = hf;
 787
 788   return 0;
 789
 790 error:
 791   if (lf) {
 792     jas_seq_destroy(lf);
 793   }
 794   if (hf) {
 795     jas_seq_destroy(hf);
 796   }
 797   return -1;
 798 }
 799
 800 #define  NNS_LIFT0(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, alpha) \
 801 { \
 802   register jpc_fix_t *lptr = (lstartptr); \
 803   register jpc_fix_t *hptr = (hstartptr); \
 804   register int n = (hendind) - (hstartind); \
 805   jpc_fix_t twoalpha = jpc_fix_mulbyint(alpha, 2); \
 806   if ((hstartind) < (lstartind)) { \
 807     jpc_fix_pluseq(*hptr, jpc_fix_mul(*lptr, (twoalpha))); \
 808     hptr += (step); \
 809     --n; \
 810   } \
 811   if ((hendind) >= (lendind)) { \
 812     --n; \
 813   } \
 814   while (n-- > 0) { \
 815     jpc_fix_pluseq(*hptr, jpc_fix_mul(jpc_fix_add(*lptr, lptr[(step)]), (alpha))); \
 816     hptr += (step); \
 817     lptr += (step); \
 818   } \
 819   if ((hendind) >= (lendind)) { \
 820     jpc_fix_pluseq(*hptr, jpc_fix_mul(*lptr, (twoalpha))); \
 821   } \
 822 }
 823
 824 #define  NNS_LIFT1(lstartptr, lstartind, lendind, hstartptr, hstartind, hendind, step, alpha) \
 825 { \
 826   register jpc_fix_t *lptr = (lstartptr); \
 827   register jpc_fix_t *hptr = (hstartptr); \
 828   register int n = (lendind) - (lstartind); \
 829   int twoalpha = jpc_fix_mulbyint(alpha, 2); \
 830   if ((hstartind) >= (lstartind)) { \
 831     jpc_fix_pluseq(*lptr, jpc_fix_mul(*hptr, (twoalpha))); \
 832     lptr += (step); \
 833     --n; \
 834   } \
 835   if ((lendind) > (hendind)) { \
 836     --n; \
 837   } \
 838   while (n-- > 0) { \
 839     jpc_fix_pluseq(*lptr, jpc_fix_mul(jpc_fix_add(*hptr, hptr[(step)]), (alpha))); \
 840     lptr += (step); \
 841     hptr += (step); \
 842   } \
 843   if ((lendind) > (hendind)) { \
 844     jpc_fix_pluseq(*lptr, jpc_fix_mul(*hptr, (twoalpha))); \
 845   } \
 846 }
 847
 848 #define  NNS_SCALE(startptr, startind, endind, step, alpha) \
 849 { \
 850   register jpc_fix_t *ptr = (startptr); \
 851   register int n = (endind) - (startind); \
 852   while (n-- > 0) { \
 853     jpc_fix_muleq(*ptr, alpha); \
 854     ptr += (step); \
 855   } \
 856 }
 857
 858 static void jpc_ns_analyze(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
 859 {
 860   jpc_fix_t *startptr;
 861   int startind;
 862   int endind;
 863   jpc_fix_t *lstartptr;
 864   int lstartind;
 865   int lendind;
 866   jpc_fix_t *hstartptr;
 867   int hstartind;
 868   int hendind;
 869   int interstep;
 870   int intrastep;
 871   int numseq;
 872
 873   /* Avoid compiler warnings about unused parameters. */
 874   qmfb = 0;
 875
 876   if (flags & JPC_QMFB1D_VERT) {
 877     interstep = 1;
 878     intrastep = jas_seq2d_rowstep(x);
 879     numseq = jas_seq2d_width(x);
 880     startind = jas_seq2d_ystart(x);
 881     endind = jas_seq2d_yend(x);
 882   } else {
 883     interstep = jas_seq2d_rowstep(x);
 884     intrastep = 1;
 885     numseq = jas_seq2d_height(x);
 886     startind = jas_seq2d_xstart(x);
 887     endind = jas_seq2d_xend(x);
 888   }
 889
 890   assert(startind < endind);
 891
 892   startptr = jas_seq2d_getref(x, jas_seq2d_xstart(x), jas_seq2d_ystart(x));
 893   if (!(flags & JPC_QMFB1D_RITIMODE)) {
 894     while (numseq-- > 0) {
 895       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 896         &lstartptr, &lstartind, &lendind, &hstartptr,
 897         &hstartind, &hendind);
 898       if (endind - startind > 1) {
 899         jpc_qmfb1d_split(startptr, startind, endind,
 900           intrastep, lstartptr, lstartind, lendind,
 901           hstartptr, hstartind, hendind);
 902         NNS_LIFT0(lstartptr, lstartind, lendind,
 903           hstartptr, hstartind, hendind, intrastep,
 904           jpc_dbltofix(-1.586134342));
 905         NNS_LIFT1(lstartptr, lstartind, lendind,
 906           hstartptr, hstartind, hendind, intrastep,
 907           jpc_dbltofix(-0.052980118));
 908         NNS_LIFT0(lstartptr, lstartind, lendind,
 909           hstartptr, hstartind, hendind, intrastep,
 910           jpc_dbltofix(0.882911075));
 911         NNS_LIFT1(lstartptr, lstartind, lendind,
 912           hstartptr, hstartind, hendind, intrastep,
 913           jpc_dbltofix(0.443506852));
 914         NNS_SCALE(lstartptr, lstartind, lendind,
 915           intrastep, jpc_dbltofix(1.0/1.23017410558578));
 916         NNS_SCALE(hstartptr, hstartind, hendind,
 917           intrastep, jpc_dbltofix(1.0/1.62578613134411));
 918       } else {
 919 #if 0
 920         if (lstartind == lendind) {
 921           *startptr = jpc_fix_asl(*startptr, 1);
 922         }
 923 #endif
 924       }
 925       startptr += interstep;
 926     }
 927   } else {
 928     /* The reversible integer-to-integer mode is not supported
 929       for this transform. */
 930     abort();
 931   }
 932 }
 933
 934 static void jpc_ns_synthesize(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
 935 {
 936   jpc_fix_t *startptr;
 937   int startind;
 938   int endind;
 939   jpc_fix_t *lstartptr;
 940   int lstartind;
 941   int lendind;
 942   jpc_fix_t *hstartptr;
 943   int hstartind;
 944   int hendind;
 945   int interstep;
 946   int intrastep;
 947   int numseq;
 948
 949   /* Avoid compiler warnings about unused parameters. */
 950   qmfb = 0;
 951
 952   if (flags & JPC_QMFB1D_VERT) {
 953     interstep = 1;
 954     intrastep = jas_seq2d_rowstep(x);
 955     numseq = jas_seq2d_width(x);
 956     startind = jas_seq2d_ystart(x);
 957     endind = jas_seq2d_yend(x);
 958   } else {
 959     interstep = jas_seq2d_rowstep(x);
 960     intrastep = 1;
 961     numseq = jas_seq2d_height(x);
 962     startind = jas_seq2d_xstart(x);
 963     endind = jas_seq2d_xend(x);
 964   }
 965
 966   assert(startind < endind);
 967
 968   startptr = jas_seq2d_getref(x, jas_seq2d_xstart(x), jas_seq2d_ystart(x));
 969   if (!(flags & JPC_QMFB1D_RITIMODE)) {
 970     while (numseq-- > 0) {
 971       jpc_qmfb1d_setup(startptr, startind, endind, intrastep,
 972         &lstartptr, &lstartind, &lendind, &hstartptr,
 973         &hstartind, &hendind);
 974       if (endind - startind > 1) {
 975         NNS_SCALE(lstartptr, lstartind, lendind,
 976           intrastep, jpc_dbltofix(1.23017410558578));
 977         NNS_SCALE(hstartptr, hstartind, hendind,
 978           intrastep, jpc_dbltofix(1.62578613134411));
 979         NNS_LIFT1(lstartptr, lstartind, lendind,
 980           hstartptr, hstartind, hendind, intrastep,
 981           jpc_dbltofix(-0.443506852));
 982         NNS_LIFT0(lstartptr, lstartind, lendind,
 983           hstartptr, hstartind, hendind, intrastep,
 984           jpc_dbltofix(-0.882911075));
 985         NNS_LIFT1(lstartptr, lstartind, lendind,
 986           hstartptr, hstartind, hendind, intrastep,
 987           jpc_dbltofix(0.052980118));
 988         NNS_LIFT0(lstartptr, lstartind, lendind,
 989           hstartptr, hstartind, hendind, intrastep,
 990           jpc_dbltofix(1.586134342));
 991         jpc_qmfb1d_join(startptr, startind, endind,
 992           intrastep, lstartptr, lstartind, lendind,
 993           hstartptr, hstartind, hendind);
 994       } else {
 995 #if 0
 996         if (lstartind == lendind) {
 997           *startptr = jpc_fix_asr(*startptr, 1);
 998         }
 999 #endif
1000       }
1001       startptr += interstep;
1002     }
1003   } else {
1004     /* The reversible integer-to-integer mode is not supported
1005       for this transform. */
1006     abort();
1007   }
1008 }
1009
1010 /******************************************************************************\
1011 *
1012 \******************************************************************************/
1013
1014 jpc_qmfb1d_t *jpc_qmfb1d_make(int qmfbid)
1015 {
1016   jpc_qmfb1d_t *qmfb;
1017   if (!(qmfb = jpc_qmfb1d_create())) {
1018     return 0;
1019   }
1020   switch (qmfbid) {
1021   case JPC_QMFB1D_FT:
1022     qmfb->ops = &jpc_ft_ops;
1023     break;
1024   case JPC_QMFB1D_NS:
1025     qmfb->ops = &jpc_ns_ops;
1026     break;
1027   default:
1028     jpc_qmfb1d_destroy(qmfb);
1029     return 0;
1030     break;
1031   }
1032   return qmfb;
1033 }
1034
1035 static jpc_qmfb1d_t *jpc_qmfb1d_create()
1036 {
1037   jpc_qmfb1d_t *qmfb;
1038   if (!(qmfb = jas_malloc(sizeof(jpc_qmfb1d_t)))) {
1039     return 0;
1040   }
1041   qmfb->ops = 0;
1042   return qmfb;
1043 }
1044
1045 jpc_qmfb1d_t *jpc_qmfb1d_copy(jpc_qmfb1d_t *qmfb)
1046 {
1047   jpc_qmfb1d_t *newqmfb;
1048
1049   if (!(newqmfb = jpc_qmfb1d_create())) {
1050     return 0;
1051   }
1052   newqmfb->ops = qmfb->ops;
1053   return newqmfb;
1054 }
1055
1056 void jpc_qmfb1d_destroy(jpc_qmfb1d_t *qmfb)
1057 {
1058   jas_free(qmfb);
1059 }
1060
1061 /******************************************************************************\
1062 *
1063 \******************************************************************************/
1064
1065 void jpc_qmfb1d_getbands(jpc_qmfb1d_t *qmfb, int flags, uint_fast32_t xstart,
1066   uint_fast32_t ystart, uint_fast32_t xend, uint_fast32_t yend, int maxbands,
1067   int *numbandsptr, jpc_qmfb1dband_t *bands)
1068 {
1069   int start;
1070   int end;
1071
1072   assert(maxbands >= 2);
1073
1074   if (flags & JPC_QMFB1D_VERT) {
1075     start = ystart;
1076     end = yend;
1077   } else {
1078     start = xstart;
1079     end = xend;
1080   }
1081   assert(jpc_qmfb1d_getnumchans(qmfb) == 2);
1082   assert(start <= end);
1083   bands[0].start = JPC_CEILDIVPOW2(start, 1);
1084   bands[0].end = JPC_CEILDIVPOW2(end, 1);
1085   bands[0].locstart = start;
1086   bands[0].locend = start + bands[0].end - bands[0].start;
1087   bands[1].start = JPC_FLOORDIVPOW2(start, 1);
1088   bands[1].end = JPC_FLOORDIVPOW2(end, 1);
1089   bands[1].locstart = bands[0].locend;
1090   bands[1].locend = bands[1].locstart + bands[1].end - bands[1].start;
1091   assert(bands[1].locend == end);
1092   *numbandsptr = 2;
1093 }
1094
1095 /******************************************************************************\
1096 *
1097 \******************************************************************************/
1098
1099 int jpc_qmfb1d_getnumchans(jpc_qmfb1d_t *qmfb)
1100 {
1101   return (*qmfb->ops->getnumchans)(qmfb);
1102 }
1103
1104 int jpc_qmfb1d_getanalfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
1105 {
1106   return (*qmfb->ops->getanalfilters)(qmfb, len, filters);
1107 }
1108
1109 int jpc_qmfb1d_getsynfilters(jpc_qmfb1d_t *qmfb, int len, jas_seq2d_t **filters)
1110 {
1111   return (*qmfb->ops->getsynfilters)(qmfb, len, filters);
1112 }
1113
1114 void jpc_qmfb1d_analyze(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
1115 {
1116   (*qmfb->ops->analyze)(qmfb, flags, x);
1117 }
1118
1119 void jpc_qmfb1d_synthesize(jpc_qmfb1d_t *qmfb, int flags, jas_seq2d_t *x)
1120 {
1121   (*qmfb->ops->synthesize)(qmfb, flags, x);
1122 }