git.maemo.org Git - opencv/blob - otherlibs/_graphics/src/libjasper/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 2007/01/15 16:09:28 vp153 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 }