Update the trunk to the OpenCV's CVS (2008-07-14)
[opencv] / 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) 2001-2006 Michael David Adams
13  * Copyright (c) 1999-2000 Image Power, Inc.
14  * Copyright (c) 1999-2000 The University of British Columbia
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.2 2008/05/26 09:40:52 vp153 Exp $
68  */
69
70 /******************************************************************************\
71 *
72 \******************************************************************************/
73
74 #undef WT_LENONE /* This is not needed due to normalization. */
75 #define WT_DOSCALE
76
77 /******************************************************************************\
78 * Includes.
79 \******************************************************************************/
80
81 #include <assert.h>
82 #include "jasper/jas_fix.h"
83 #include "jasper/jas_malloc.h"
84 #include "jasper/jas_math.h"
85
86 #include "jpc_qmfb.h"
87 #include "jpc_tsfb.h"
88 #include "jpc_math.h"
89
90 /******************************************************************************\
91 *
92 \******************************************************************************/
93
94 #define QMFB_SPLITBUFSIZE       4096
95 #define QMFB_JOINBUFSIZE        4096
96
97 int jpc_ft_analyze(jpc_fix_t *a, int xstart, int ystart, int width, int height,
98   int stride);
99 int jpc_ft_synthesize(int *a, int xstart, int ystart, int width, int height,
100   int stride);
101
102 int jpc_ns_analyze(jpc_fix_t *a, int xstart, int ystart, int width, int height,
103   int stride);
104 int jpc_ns_synthesize(jpc_fix_t *a, int xstart, int ystart, int width,
105   int height, int stride);
106
107 void jpc_ft_fwdlift_row(jpc_fix_t *a, int numcols, int parity);
108 void jpc_ft_fwdlift_col(jpc_fix_t *a, int numrows, int stride,
109   int parity);
110 void jpc_ft_fwdlift_colgrp(jpc_fix_t *a, int numrows, int stride,
111   int parity);
112 void jpc_ft_fwdlift_colres(jpc_fix_t *a, int numrows, int numcols,
113   int stride, int parity);
114
115 void jpc_ft_invlift_row(jpc_fix_t *a, int numcols, int parity);
116 void jpc_ft_invlift_col(jpc_fix_t *a, int numrows, int stride,
117   int parity);
118 void jpc_ft_invlift_colgrp(jpc_fix_t *a, int numrows, int stride,
119   int parity);
120 void jpc_ft_invlift_colres(jpc_fix_t *a, int numrows, int numcols,
121   int stride, int parity);
122
123 void jpc_ns_fwdlift_row(jpc_fix_t *a, int numcols, int parity);
124 void jpc_ns_fwdlift_colgrp(jpc_fix_t *a, int numrows, int stride, int parity);
125 void jpc_ns_fwdlift_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
126   int parity);
127 void jpc_ns_invlift_row(jpc_fix_t *a, int numcols, int parity);
128 void jpc_ns_invlift_colgrp(jpc_fix_t *a, int numrows, int stride, int parity);
129 void jpc_ns_invlift_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
130   int parity);
131
132 void jpc_qmfb_split_row(jpc_fix_t *a, int numcols, int parity);
133 void jpc_qmfb_split_col(jpc_fix_t *a, int numrows, int stride, int parity);
134 void jpc_qmfb_split_colgrp(jpc_fix_t *a, int numrows, int stride, int parity);
135 void jpc_qmfb_split_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
136   int parity);
137
138 void jpc_qmfb_join_row(jpc_fix_t *a, int numcols, int parity);
139 void jpc_qmfb_join_col(jpc_fix_t *a, int numrows, int stride, int parity);
140 void jpc_qmfb_join_colgrp(jpc_fix_t *a, int numrows, int stride, int parity);
141 void jpc_qmfb_join_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
142   int parity);
143
144 double jpc_ft_lpenergywts[32] = {
145         1.2247448713915889,
146         1.6583123951776999,
147         2.3184046238739260,
148         3.2691742076555053,
149         4.6199296531440819,
150         6.5323713152269596,
151         9.2377452606141937,
152         13.0639951297449581,
153         18.4752262333915667,
154         26.1278968190610392,
155         36.9504194305524791,
156         52.2557819580462777,
157         73.9008347315741645,
158         104.5115624560829133,
159         147.8016689469569656,
160         209.0231247296646018,
161         295.6033378293900000,
162         418.0462494347059419,
163         591.2066756503630813,
164         836.0924988714708661,
165         /* approximations */
166         836.0924988714708661,
167         836.0924988714708661,
168         836.0924988714708661,
169         836.0924988714708661,
170         836.0924988714708661,
171         836.0924988714708661,
172         836.0924988714708661,
173         836.0924988714708661,
174         836.0924988714708661,
175         836.0924988714708661,
176         836.0924988714708661,
177         836.0924988714708661
178 };
179
180 double jpc_ft_hpenergywts[32] = {
181         0.8477912478906585,
182         0.9601432184835760,
183         1.2593401049756179,
184         1.7444107171191079,
185         2.4538713036750726,
186         3.4656517695088755,
187         4.8995276398597856,
188         6.9283970402160842,
189         9.7980274940131444,
190         13.8564306871112652,
191         19.5959265076535587,
192         27.7128159494245487,
193         39.1918369552045860,
194         55.4256262207444053,
195         78.3836719028959124,
196         110.8512517317256822,
197         156.7673435548526868,
198         221.7025033739244293,
199         313.5346870787551552,
200         443.4050067351659550,
201         /* approximations */
202         443.4050067351659550,
203         443.4050067351659550,
204         443.4050067351659550,
205         443.4050067351659550,
206         443.4050067351659550,
207         443.4050067351659550,
208         443.4050067351659550,
209         443.4050067351659550,
210         443.4050067351659550,
211         443.4050067351659550,
212         443.4050067351659550,
213         443.4050067351659550
214 };
215
216 double jpc_ns_lpenergywts[32] = {
217         1.4021081679297411,
218         2.0303718560817923,
219         2.9011625562785555,
220         4.1152851751758002,
221         5.8245108637728071,
222         8.2387599345725171,
223         11.6519546479210838,
224         16.4785606470644375,
225         23.3042776444606794,
226         32.9572515613740435,
227         46.6086013487782793,
228         65.9145194076860861,
229         93.2172084551803977,
230         131.8290408510004283,
231         186.4344176300625691,
232         263.6580819564562148,
233         372.8688353500955373,
234         527.3161639447193920,
235         745.7376707114038936,
236         1054.6323278917823245,
237         /* approximations follow */
238         1054.6323278917823245,
239         1054.6323278917823245,
240         1054.6323278917823245,
241         1054.6323278917823245,
242         1054.6323278917823245,
243         1054.6323278917823245,
244         1054.6323278917823245,
245         1054.6323278917823245,
246         1054.6323278917823245,
247         1054.6323278917823245,
248         1054.6323278917823245,
249         1054.6323278917823245
250 };
251
252 double jpc_ns_hpenergywts[32] = {
253         1.4425227650161456,
254         1.9669426082455688,
255         2.8839248082788891,
256         4.1475208393432981,
257         5.8946497530677817,
258         8.3471789178590949,
259         11.8086046551047463,
260         16.7012780415647804,
261         23.6196657032246620,
262         33.4034255108592362,
263         47.2396388881632632,
264         66.8069597416714061,
265         94.4793162154500692,
266         133.6139330736999113,
267         188.9586372358249378,
268         267.2278678461869390,
269         377.9172750722391356,
270         534.4557359047058753,
271         755.8345502191498326,
272         1068.9114718353569060,
273         /* approximations follow */
274         1068.9114718353569060,
275         1068.9114718353569060,
276         1068.9114718353569060,
277         1068.9114718353569060,
278         1068.9114718353569060,
279         1068.9114718353569060,
280         1068.9114718353569060,
281         1068.9114718353569060,
282         1068.9114718353569060,
283         1068.9114718353569060,
284         1068.9114718353569060
285 };
286
287 jpc_qmfb2d_t jpc_ft_qmfb2d = {
288         jpc_ft_analyze,
289         jpc_ft_synthesize,
290         jpc_ft_lpenergywts,
291         jpc_ft_hpenergywts
292 };
293
294 jpc_qmfb2d_t jpc_ns_qmfb2d = {
295         jpc_ns_analyze,
296         jpc_ns_synthesize,
297         jpc_ns_lpenergywts,
298         jpc_ns_hpenergywts
299 };
300
301 /******************************************************************************\
302 * generic
303 \******************************************************************************/
304
305 void jpc_qmfb_split_row(jpc_fix_t *a, int numcols, int parity)
306 {
307
308         int bufsize = JPC_CEILDIVPOW2(numcols, 1);
309 #if !defined(HAVE_VLA)
310         jpc_fix_t splitbuf[QMFB_SPLITBUFSIZE];
311 #else
312         jpc_fix_t splitbuf[bufsize];
313 #endif
314         jpc_fix_t *buf = splitbuf;
315         register jpc_fix_t *srcptr;
316         register jpc_fix_t *dstptr;
317         register int n;
318         register int m;
319         int hstartcol;
320
321 #if !defined(HAVE_VLA)
322         /* Get a buffer. */
323         if (bufsize > QMFB_SPLITBUFSIZE) {
324                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
325                         /* We have no choice but to commit suicide in this case. */
326                         abort();
327                 }
328         }
329 #endif
330
331         if (numcols >= 2) {
332                 hstartcol = (numcols + 1 - parity) >> 1;
333                 m = (parity) ? hstartcol : (numcols - hstartcol);
334                 /* Save the samples destined for the highpass channel. */
335                 n = m;
336                 dstptr = buf;
337                 srcptr = &a[1 - parity];
338                 while (n-- > 0) {
339                         *dstptr = *srcptr;
340                         ++dstptr;
341                         srcptr += 2;
342                 }
343                 /* Copy the appropriate samples into the lowpass channel. */
344                 dstptr = &a[1 - parity];
345                 srcptr = &a[2 - parity];
346                 n = numcols - m - (!parity);
347                 while (n-- > 0) {
348                         *dstptr = *srcptr;
349                         ++dstptr;
350                         srcptr += 2;
351                 }
352                 /* Copy the saved samples into the highpass channel. */
353                 dstptr = &a[hstartcol];
354                 srcptr = buf;
355                 n = m;
356                 while (n-- > 0) {
357                         *dstptr = *srcptr;
358                         ++dstptr;
359                         ++srcptr;
360                 }
361         }
362
363 #if !defined(HAVE_VLA)
364         /* If the split buffer was allocated on the heap, free this memory. */
365         if (buf != splitbuf) {
366                 jas_free(buf);
367         }
368 #endif
369
370 }
371
372 void jpc_qmfb_split_col(jpc_fix_t *a, int numrows, int stride,
373   int parity)
374 {
375
376         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
377 #if !defined(HAVE_VLA)
378         jpc_fix_t splitbuf[QMFB_SPLITBUFSIZE];
379 #else
380         jpc_fix_t splitbuf[bufsize];
381 #endif
382         jpc_fix_t *buf = splitbuf;
383         register jpc_fix_t *srcptr;
384         register jpc_fix_t *dstptr;
385         register int n;
386         register int m;
387         int hstartcol;
388
389 #if !defined(HAVE_VLA)
390         /* Get a buffer. */
391         if (bufsize > QMFB_SPLITBUFSIZE) {
392                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
393                         /* We have no choice but to commit suicide in this case. */
394                         abort();
395                 }
396         }
397 #endif
398
399         if (numrows >= 2) {
400                 hstartcol = (numrows + 1 - parity) >> 1;
401                 m = (parity) ? hstartcol : (numrows - hstartcol);
402                 /* Save the samples destined for the highpass channel. */
403                 n = m;
404                 dstptr = buf;
405                 srcptr = &a[(1 - parity) * stride];
406                 while (n-- > 0) {
407                         *dstptr = *srcptr;
408                         ++dstptr;
409                         srcptr += stride << 1;
410                 }
411                 /* Copy the appropriate samples into the lowpass channel. */
412                 dstptr = &a[(1 - parity) * stride];
413                 srcptr = &a[(2 - parity) * stride];
414                 n = numrows - m - (!parity);
415                 while (n-- > 0) {
416                         *dstptr = *srcptr;
417                         dstptr += stride;
418                         srcptr += stride << 1;
419                 }
420                 /* Copy the saved samples into the highpass channel. */
421                 dstptr = &a[hstartcol * stride];
422                 srcptr = buf;
423                 n = m;
424                 while (n-- > 0) {
425                         *dstptr = *srcptr;
426                         dstptr += stride;
427                         ++srcptr;
428                 }
429         }
430
431 #if !defined(HAVE_VLA)
432         /* If the split buffer was allocated on the heap, free this memory. */
433         if (buf != splitbuf) {
434                 jas_free(buf);
435         }
436 #endif
437
438 }
439
440 void jpc_qmfb_split_colgrp(jpc_fix_t *a, int numrows, int stride,
441   int parity)
442 {
443
444         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
445 #if !defined(HAVE_VLA)
446         jpc_fix_t splitbuf[QMFB_SPLITBUFSIZE * JPC_QMFB_COLGRPSIZE];
447 #else
448         jpc_fix_t splitbuf[bufsize * JPC_QMFB_COLGRPSIZE];
449 #endif
450         jpc_fix_t *buf = splitbuf;
451         jpc_fix_t *srcptr;
452         jpc_fix_t *dstptr;
453         register jpc_fix_t *srcptr2;
454         register jpc_fix_t *dstptr2;
455         register int n;
456         register int i;
457         int m;
458         int hstartcol;
459
460 #if !defined(HAVE_VLA)
461         /* Get a buffer. */
462         if (bufsize > QMFB_SPLITBUFSIZE) {
463                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
464                         /* We have no choice but to commit suicide in this case. */
465                         abort();
466                 }
467         }
468 #endif
469
470         if (numrows >= 2) {
471                 hstartcol = (numrows + 1 - parity) >> 1;
472                 m = (parity) ? hstartcol : (numrows - hstartcol);
473                 /* Save the samples destined for the highpass channel. */
474                 n = m;
475                 dstptr = buf;
476                 srcptr = &a[(1 - parity) * stride];
477                 while (n-- > 0) {
478                         dstptr2 = dstptr;
479                         srcptr2 = srcptr;
480                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
481                                 *dstptr2 = *srcptr2;
482                                 ++dstptr2;
483                                 ++srcptr2;
484                         }
485                         dstptr += JPC_QMFB_COLGRPSIZE;
486                         srcptr += stride << 1;
487                 }
488                 /* Copy the appropriate samples into the lowpass channel. */
489                 dstptr = &a[(1 - parity) * stride];
490                 srcptr = &a[(2 - parity) * stride];
491                 n = numrows - m - (!parity);
492                 while (n-- > 0) {
493                         dstptr2 = dstptr;
494                         srcptr2 = srcptr;
495                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
496                                 *dstptr2 = *srcptr2;
497                                 ++dstptr2;
498                                 ++srcptr2;
499                         }
500                         dstptr += stride;
501                         srcptr += stride << 1;
502                 }
503                 /* Copy the saved samples into the highpass channel. */
504                 dstptr = &a[hstartcol * stride];
505                 srcptr = buf;
506                 n = m;
507                 while (n-- > 0) {
508                         dstptr2 = dstptr;
509                         srcptr2 = srcptr;
510                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
511                                 *dstptr2 = *srcptr2;
512                                 ++dstptr2;
513                                 ++srcptr2;
514                         }
515                         dstptr += stride;
516                         srcptr += JPC_QMFB_COLGRPSIZE;
517                 }
518         }
519
520 #if !defined(HAVE_VLA)
521         /* If the split buffer was allocated on the heap, free this memory. */
522         if (buf != splitbuf) {
523                 jas_free(buf);
524         }
525 #endif
526
527 }
528
529 void jpc_qmfb_split_colres(jpc_fix_t *a, int numrows, int numcols,
530   int stride, int parity)
531 {
532
533         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
534 #if !defined(HAVE_VLA)
535         jpc_fix_t splitbuf[QMFB_SPLITBUFSIZE * JPC_QMFB_COLGRPSIZE];
536 #else
537         jpc_fix_t splitbuf[bufsize * numcols];
538 #endif
539         jpc_fix_t *buf = splitbuf;
540         jpc_fix_t *srcptr;
541         jpc_fix_t *dstptr;
542         register jpc_fix_t *srcptr2;
543         register jpc_fix_t *dstptr2;
544         register int n;
545         register int i;
546         int m;
547         int hstartcol;
548
549 #if !defined(HAVE_VLA)
550         /* Get a buffer. */
551         if (bufsize > QMFB_SPLITBUFSIZE) {
552                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
553                         /* We have no choice but to commit suicide in this case. */
554                         abort();
555                 }
556         }
557 #endif
558
559         if (numrows >= 2) {
560                 hstartcol = (numrows + 1 - parity) >> 1;
561                 m = (parity) ? hstartcol : (numrows - hstartcol);
562                 /* Save the samples destined for the highpass channel. */
563                 n = m;
564                 dstptr = buf;
565                 srcptr = &a[(1 - parity) * stride];
566                 while (n-- > 0) {
567                         dstptr2 = dstptr;
568                         srcptr2 = srcptr;
569                         for (i = 0; i < numcols; ++i) {
570                                 *dstptr2 = *srcptr2;
571                                 ++dstptr2;
572                                 ++srcptr2;
573                         }
574                         dstptr += numcols;
575                         srcptr += stride << 1;
576                 }
577                 /* Copy the appropriate samples into the lowpass channel. */
578                 dstptr = &a[(1 - parity) * stride];
579                 srcptr = &a[(2 - parity) * stride];
580                 n = numrows - m - (!parity);
581                 while (n-- > 0) {
582                         dstptr2 = dstptr;
583                         srcptr2 = srcptr;
584                         for (i = 0; i < numcols; ++i) {
585                                 *dstptr2 = *srcptr2;
586                                 ++dstptr2;
587                                 ++srcptr2;
588                         }
589                         dstptr += stride;
590                         srcptr += stride << 1;
591                 }
592                 /* Copy the saved samples into the highpass channel. */
593                 dstptr = &a[hstartcol * stride];
594                 srcptr = buf;
595                 n = m;
596                 while (n-- > 0) {
597                         dstptr2 = dstptr;
598                         srcptr2 = srcptr;
599                         for (i = 0; i < numcols; ++i) {
600                                 *dstptr2 = *srcptr2;
601                                 ++dstptr2;
602                                 ++srcptr2;
603                         }
604                         dstptr += stride;
605                         srcptr += numcols;
606                 }
607         }
608
609 #if !defined(HAVE_VLA)
610         /* If the split buffer was allocated on the heap, free this memory. */
611         if (buf != splitbuf) {
612                 jas_free(buf);
613         }
614 #endif
615
616 }
617
618 void jpc_qmfb_join_row(jpc_fix_t *a, int numcols, int parity)
619 {
620
621         int bufsize = JPC_CEILDIVPOW2(numcols, 1);
622 #if !defined(HAVE_VLA)
623         jpc_fix_t joinbuf[QMFB_JOINBUFSIZE];
624 #else
625         jpc_fix_t joinbuf[bufsize];
626 #endif
627         jpc_fix_t *buf = joinbuf;
628         register jpc_fix_t *srcptr;
629         register jpc_fix_t *dstptr;
630         register int n;
631         int hstartcol;
632
633 #if !defined(HAVE_VLA)
634         /* Allocate memory for the join buffer from the heap. */
635         if (bufsize > QMFB_JOINBUFSIZE) {
636                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
637                         /* We have no choice but to commit suicide. */
638                         abort();
639                 }
640         }
641 #endif
642
643         hstartcol = (numcols + 1 - parity) >> 1;
644
645         /* Save the samples from the lowpass channel. */
646         n = hstartcol;
647         srcptr = &a[0];
648         dstptr = buf;
649         while (n-- > 0) {
650                 *dstptr = *srcptr;
651                 ++srcptr;
652                 ++dstptr;
653         }
654         /* Copy the samples from the highpass channel into place. */
655         srcptr = &a[hstartcol];
656         dstptr = &a[1 - parity];
657         n = numcols - hstartcol;
658         while (n-- > 0) {
659                 *dstptr = *srcptr;
660                 dstptr += 2;
661                 ++srcptr;
662         }
663         /* Copy the samples from the lowpass channel into place. */
664         srcptr = buf;
665         dstptr = &a[parity];
666         n = hstartcol;
667         while (n-- > 0) {
668                 *dstptr = *srcptr;
669                 dstptr += 2;
670                 ++srcptr;
671         }
672
673 #if !defined(HAVE_VLA)
674         /* If the join buffer was allocated on the heap, free this memory. */
675         if (buf != joinbuf) {
676                 jas_free(buf);
677         }
678 #endif
679
680 }
681
682 void jpc_qmfb_join_col(jpc_fix_t *a, int numrows, int stride,
683   int parity)
684 {
685
686         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
687 #if !defined(HAVE_VLA)
688         jpc_fix_t joinbuf[QMFB_JOINBUFSIZE];
689 #else
690         jpc_fix_t joinbuf[bufsize];
691 #endif
692         jpc_fix_t *buf = joinbuf;
693         register jpc_fix_t *srcptr;
694         register jpc_fix_t *dstptr;
695         register int n;
696         int hstartcol;
697
698 #if !defined(HAVE_VLA)
699         /* Allocate memory for the join buffer from the heap. */
700         if (bufsize > QMFB_JOINBUFSIZE) {
701                 if (!(buf = jas_malloc(bufsize * sizeof(jpc_fix_t)))) {
702                         /* We have no choice but to commit suicide. */
703                         abort();
704                 }
705         }
706 #endif
707
708         hstartcol = (numrows + 1 - parity) >> 1;
709
710         /* Save the samples from the lowpass channel. */
711         n = hstartcol;
712         srcptr = &a[0];
713         dstptr = buf;
714         while (n-- > 0) {
715                 *dstptr = *srcptr;
716                 srcptr += stride;
717                 ++dstptr;
718         }
719         /* Copy the samples from the highpass channel into place. */
720         srcptr = &a[hstartcol * stride];
721         dstptr = &a[(1 - parity) * stride];
722         n = numrows - hstartcol;
723         while (n-- > 0) {
724                 *dstptr = *srcptr;
725                 dstptr += 2 * stride;
726                 srcptr += stride;
727         }
728         /* Copy the samples from the lowpass channel into place. */
729         srcptr = buf;
730         dstptr = &a[parity * stride];
731         n = hstartcol;
732         while (n-- > 0) {
733                 *dstptr = *srcptr;
734                 dstptr += 2 * stride;
735                 ++srcptr;
736         }
737
738 #if !defined(HAVE_VLA)
739         /* If the join buffer was allocated on the heap, free this memory. */
740         if (buf != joinbuf) {
741                 jas_free(buf);
742         }
743 #endif
744
745 }
746
747 void jpc_qmfb_join_colgrp(jpc_fix_t *a, int numrows, int stride,
748   int parity)
749 {
750
751         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
752 #if !defined(HAVE_VLA)
753         jpc_fix_t joinbuf[QMFB_JOINBUFSIZE * JPC_QMFB_COLGRPSIZE];
754 #else
755         jpc_fix_t joinbuf[bufsize * JPC_QMFB_COLGRPSIZE];
756 #endif
757         jpc_fix_t *buf = joinbuf;
758         jpc_fix_t *srcptr;
759         jpc_fix_t *dstptr;
760         register jpc_fix_t *srcptr2;
761         register jpc_fix_t *dstptr2;
762         register int n;
763         register int i;
764         int hstartcol;
765
766 #if !defined(HAVE_VLA)
767         /* Allocate memory for the join buffer from the heap. */
768         if (bufsize > QMFB_JOINBUFSIZE) {
769                 if (!(buf = jas_malloc(bufsize * JPC_QMFB_COLGRPSIZE * sizeof(jpc_fix_t)))) {
770                         /* We have no choice but to commit suicide. */
771                         abort();
772                 }
773         }
774 #endif
775
776         hstartcol = (numrows + 1 - parity) >> 1;
777
778         /* Save the samples from the lowpass channel. */
779         n = hstartcol;
780         srcptr = &a[0];
781         dstptr = buf;
782         while (n-- > 0) {
783                 dstptr2 = dstptr;
784                 srcptr2 = srcptr;
785                 for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
786                         *dstptr2 = *srcptr2;
787                         ++dstptr2;
788                         ++srcptr2;
789                 }
790                 srcptr += stride;
791                 dstptr += JPC_QMFB_COLGRPSIZE;
792         }
793         /* Copy the samples from the highpass channel into place. */
794         srcptr = &a[hstartcol * stride];
795         dstptr = &a[(1 - parity) * stride];
796         n = numrows - hstartcol;
797         while (n-- > 0) {
798                 dstptr2 = dstptr;
799                 srcptr2 = srcptr;
800                 for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
801                         *dstptr2 = *srcptr2;
802                         ++dstptr2;
803                         ++srcptr2;
804                 }
805                 dstptr += 2 * stride;
806                 srcptr += stride;
807         }
808         /* Copy the samples from the lowpass channel into place. */
809         srcptr = buf;
810         dstptr = &a[parity * stride];
811         n = hstartcol;
812         while (n-- > 0) {
813                 dstptr2 = dstptr;
814                 srcptr2 = srcptr;
815                 for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
816                         *dstptr2 = *srcptr2;
817                         ++dstptr2;
818                         ++srcptr2;
819                 }
820                 dstptr += 2 * stride;
821                 srcptr += JPC_QMFB_COLGRPSIZE;
822         }
823
824 #if !defined(HAVE_VLA)
825         /* If the join buffer was allocated on the heap, free this memory. */
826         if (buf != joinbuf) {
827                 jas_free(buf);
828         }
829 #endif
830
831 }
832
833 void jpc_qmfb_join_colres(jpc_fix_t *a, int numrows, int numcols,
834   int stride, int parity)
835 {
836
837         int bufsize = JPC_CEILDIVPOW2(numrows, 1);
838 #if !defined(HAVE_VLA)
839         jpc_fix_t joinbuf[QMFB_JOINBUFSIZE * JPC_QMFB_COLGRPSIZE];
840 #else
841         jpc_fix_t joinbuf[bufsize * numcols];
842 #endif
843         jpc_fix_t *buf = joinbuf;
844         jpc_fix_t *srcptr;
845         jpc_fix_t *dstptr;
846         register jpc_fix_t *srcptr2;
847         register jpc_fix_t *dstptr2;
848         register int n;
849         register int i;
850         int hstartcol;
851
852 #if !defined(HAVE_VLA)
853         /* Allocate memory for the join buffer from the heap. */
854         if (bufsize > QMFB_JOINBUFSIZE) {
855                 if (!(buf = jas_malloc(bufsize * numcols * sizeof(jpc_fix_t)))) {
856                         /* We have no choice but to commit suicide. */
857                         abort();
858                 }
859         }
860 #endif
861
862         hstartcol = (numrows + 1 - parity) >> 1;
863
864         /* Save the samples from the lowpass channel. */
865         n = hstartcol;
866         srcptr = &a[0];
867         dstptr = buf;
868         while (n-- > 0) {
869                 dstptr2 = dstptr;
870                 srcptr2 = srcptr;
871                 for (i = 0; i < numcols; ++i) {
872                         *dstptr2 = *srcptr2;
873                         ++dstptr2;
874                         ++srcptr2;
875                 }
876                 srcptr += stride;
877                 dstptr += numcols;
878         }
879         /* Copy the samples from the highpass channel into place. */
880         srcptr = &a[hstartcol * stride];
881         dstptr = &a[(1 - parity) * stride];
882         n = numrows - hstartcol;
883         while (n-- > 0) {
884                 dstptr2 = dstptr;
885                 srcptr2 = srcptr;
886                 for (i = 0; i < numcols; ++i) {
887                         *dstptr2 = *srcptr2;
888                         ++dstptr2;
889                         ++srcptr2;
890                 }
891                 dstptr += 2 * stride;
892                 srcptr += stride;
893         }
894         /* Copy the samples from the lowpass channel into place. */
895         srcptr = buf;
896         dstptr = &a[parity * stride];
897         n = hstartcol;
898         while (n-- > 0) {
899                 dstptr2 = dstptr;
900                 srcptr2 = srcptr;
901                 for (i = 0; i < numcols; ++i) {
902                         *dstptr2 = *srcptr2;
903                         ++dstptr2;
904                         ++srcptr2;
905                 }
906                 dstptr += 2 * stride;
907                 srcptr += numcols;
908         }
909
910 #if !defined(HAVE_VLA)
911         /* If the join buffer was allocated on the heap, free this memory. */
912         if (buf != joinbuf) {
913                 jas_free(buf);
914         }
915 #endif
916
917 }
918
919 /******************************************************************************\
920 * 5/3 transform
921 \******************************************************************************/
922
923 void jpc_ft_fwdlift_row(jpc_fix_t *a, int numcols, int parity)
924 {
925
926         register jpc_fix_t *lptr;
927         register jpc_fix_t *hptr;
928         register int n;
929         int llen;
930
931         llen = (numcols + 1 - parity) >> 1;
932
933         if (numcols > 1) {
934
935                 /* Apply the first lifting step. */
936                 lptr = &a[0];
937                 hptr = &a[llen];
938                 if (parity) {
939                         hptr[0] -= lptr[0];
940                         ++hptr;
941                 }
942                 n = numcols - llen - parity - (parity == (numcols & 1));
943                 while (n-- > 0) {
944                         hptr[0] -= (lptr[0] + lptr[1]) >> 1;
945                         ++hptr;
946                         ++lptr;
947                 }
948                 if (parity == (numcols & 1)) {
949                         hptr[0] -= lptr[0];
950                 }
951
952                 /* Apply the second lifting step. */
953                 lptr = &a[0];
954                 hptr = &a[llen];
955                 if (!parity) {
956                         lptr[0] += (hptr[0] + 1) >> 1;
957                         ++lptr;
958                 }
959                 n = llen - (!parity) - (parity != (numcols & 1));
960                 while (n-- > 0) {
961                         lptr[0] += (hptr[0] + hptr[1] + 2) >> 2;
962                         ++lptr;
963                         ++hptr;
964                 }
965                 if (parity != (numcols & 1)) {
966                         lptr[0] += (hptr[0] + 1) >> 1;
967                 }
968
969         } else {
970
971                 if (parity) {
972                         lptr = &a[0];
973                         lptr[0] <<= 1;
974                 }
975
976         }
977
978 }
979
980 void jpc_ft_fwdlift_col(jpc_fix_t *a, int numrows, int stride, int parity)
981 {
982
983         jpc_fix_t *lptr;
984         jpc_fix_t *hptr;
985 #if 0
986         register jpc_fix_t *lptr2;
987         register jpc_fix_t *hptr2;
988         register int i;
989 #endif
990         register int n;
991         int llen;
992
993         llen = (numrows + 1 - parity) >> 1;
994
995         if (numrows > 1) {
996
997                 /* Apply the first lifting step. */
998                 lptr = &a[0];
999                 hptr = &a[llen * stride];
1000                 if (parity) {
1001                         hptr[0] -= lptr[0];
1002                         hptr += stride;
1003                 }
1004                 n = numrows - llen - parity - (parity == (numrows & 1));
1005                 while (n-- > 0) {
1006                         hptr[0] -= (lptr[0] + lptr[stride]) >> 1;
1007                         hptr += stride;
1008                         lptr += stride;
1009                 }
1010                 if (parity == (numrows & 1)) {
1011                         hptr[0] -= lptr[0];
1012                 }
1013
1014                 /* Apply the second lifting step. */
1015                 lptr = &a[0];
1016                 hptr = &a[llen * stride];
1017                 if (!parity) {
1018                         lptr[0] += (hptr[0] + 1) >> 1;
1019                         lptr += stride;
1020                 }
1021                 n = llen - (!parity) - (parity != (numrows & 1));
1022                 while (n-- > 0) {
1023                         lptr[0] += (hptr[0] + hptr[stride] + 2) >> 2;
1024                         lptr += stride;
1025                         hptr += stride;
1026                 }
1027                 if (parity != (numrows & 1)) {
1028                         lptr[0] += (hptr[0] + 1) >> 1;
1029                 }
1030
1031         } else {
1032
1033                 if (parity) {
1034                         lptr = &a[0];
1035                         lptr[0] <<= 1;
1036                 }
1037
1038         }
1039
1040 }
1041
1042 void jpc_ft_fwdlift_colgrp(jpc_fix_t *a, int numrows, int stride, int parity)
1043 {
1044
1045         jpc_fix_t *lptr;
1046         jpc_fix_t *hptr;
1047         register jpc_fix_t *lptr2;
1048         register jpc_fix_t *hptr2;
1049         register int n;
1050         register int i;
1051         int llen;
1052
1053         llen = (numrows + 1 - parity) >> 1;
1054
1055         if (numrows > 1) {
1056
1057                 /* Apply the first lifting step. */
1058                 lptr = &a[0];
1059                 hptr = &a[llen * stride];
1060                 if (parity) {
1061                         lptr2 = lptr;
1062                         hptr2 = hptr;
1063                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1064                                 hptr2[0] -= lptr2[0];
1065                                 ++hptr2;
1066                                 ++lptr2;
1067                         }
1068                         hptr += stride;
1069                 }
1070                 n = numrows - llen - parity - (parity == (numrows & 1));
1071                 while (n-- > 0) {
1072                         lptr2 = lptr;
1073                         hptr2 = hptr;
1074                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1075                                 hptr2[0] -= (lptr2[0] + lptr2[stride]) >> 1;
1076                                 ++lptr2;
1077                                 ++hptr2;
1078                         }
1079                         hptr += stride;
1080                         lptr += stride;
1081                 }
1082                 if (parity == (numrows & 1)) {
1083                         lptr2 = lptr;
1084                         hptr2 = hptr;
1085                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1086                                 hptr2[0] -= lptr2[0];
1087                                 ++lptr2;
1088                                 ++hptr2;
1089                         }
1090                 }
1091
1092                 /* Apply the second lifting step. */
1093                 lptr = &a[0];
1094                 hptr = &a[llen * stride];
1095                 if (!parity) {
1096                         lptr2 = lptr;
1097                         hptr2 = hptr;
1098                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1099                                 lptr2[0] += (hptr2[0] + 1) >> 1;
1100                                 ++lptr2;
1101                                 ++hptr2;
1102                         }
1103                         lptr += stride;
1104                 }
1105                 n = llen - (!parity) - (parity != (numrows & 1));
1106                 while (n-- > 0) {
1107                         lptr2 = lptr;
1108                         hptr2 = hptr;
1109                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1110                                 lptr2[0] += (hptr2[0] + hptr2[stride] + 2) >> 2;
1111                                 ++lptr2;
1112                                 ++hptr2;
1113                         }
1114                         lptr += stride;
1115                         hptr += stride;
1116                 }
1117                 if (parity != (numrows & 1)) {
1118                         lptr2 = lptr;
1119                         hptr2 = hptr;
1120                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1121                                 lptr2[0] += (hptr2[0] + 1) >> 1;
1122                                 ++lptr2;
1123                                 ++hptr2;
1124                         }
1125                 }
1126
1127         } else {
1128
1129                 if (parity) {
1130                         lptr2 = &a[0];
1131                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1132                                 lptr2[0] <<= 1;
1133                                 ++lptr2;
1134                         }
1135                 }
1136
1137         }
1138
1139 }
1140
1141 void jpc_ft_fwdlift_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
1142   int parity)
1143 {
1144
1145         jpc_fix_t *lptr;
1146         jpc_fix_t *hptr;
1147         register jpc_fix_t *lptr2;
1148         register jpc_fix_t *hptr2;
1149         register int n;
1150         register int i;
1151         int llen;
1152
1153         llen = (numrows + 1 - parity) >> 1;
1154
1155         if (numrows > 1) {
1156
1157                 /* Apply the first lifting step. */
1158                 lptr = &a[0];
1159                 hptr = &a[llen * stride];
1160                 if (parity) {
1161                         lptr2 = lptr;
1162                         hptr2 = hptr;
1163                         for (i = 0; i < numcols; ++i) {
1164                                 hptr2[0] -= lptr2[0];
1165                                 ++hptr2;
1166                                 ++lptr2;
1167                         }
1168                         hptr += stride;
1169                 }
1170                 n = numrows - llen - parity - (parity == (numrows & 1));
1171                 while (n-- > 0) {
1172                         lptr2 = lptr;
1173                         hptr2 = hptr;
1174                         for (i = 0; i < numcols; ++i) {
1175                                 hptr2[0] -= (lptr2[0] + lptr2[stride]) >> 1;
1176                                 ++lptr2;
1177                                 ++hptr2;
1178                         }
1179                         hptr += stride;
1180                         lptr += stride;
1181                 }
1182                 if (parity == (numrows & 1)) {
1183                         lptr2 = lptr;
1184                         hptr2 = hptr;
1185                         for (i = 0; i < numcols; ++i) {
1186                                 hptr2[0] -= lptr2[0];
1187                                 ++lptr2;
1188                                 ++hptr2;
1189                         }
1190                 }
1191
1192                 /* Apply the second lifting step. */
1193                 lptr = &a[0];
1194                 hptr = &a[llen * stride];
1195                 if (!parity) {
1196                         lptr2 = lptr;
1197                         hptr2 = hptr;
1198                         for (i = 0; i < numcols; ++i) {
1199                                 lptr2[0] += (hptr2[0] + 1) >> 1;
1200                                 ++lptr2;
1201                                 ++hptr2;
1202                         }
1203                         lptr += stride;
1204                 }
1205                 n = llen - (!parity) - (parity != (numrows & 1));
1206                 while (n-- > 0) {
1207                         lptr2 = lptr;
1208                         hptr2 = hptr;
1209                         for (i = 0; i < numcols; ++i) {
1210                                 lptr2[0] += (hptr2[0] + hptr2[stride] + 2) >> 2;
1211                                 ++lptr2;
1212                                 ++hptr2;
1213                         }
1214                         lptr += stride;
1215                         hptr += stride;
1216                 }
1217                 if (parity != (numrows & 1)) {
1218                         lptr2 = lptr;
1219                         hptr2 = hptr;
1220                         for (i = 0; i < numcols; ++i) {
1221                                 lptr2[0] += (hptr2[0] + 1) >> 1;
1222                                 ++lptr2;
1223                                 ++hptr2;
1224                         }
1225                 }
1226
1227         } else {
1228
1229                 if (parity) {
1230                         lptr2 = &a[0];
1231                         for (i = 0; i < numcols; ++i) {
1232                                 lptr2[0] <<= 1;
1233                                 ++lptr2;
1234                         }
1235                 }
1236
1237         }
1238
1239 }
1240
1241 void jpc_ft_invlift_row(jpc_fix_t *a, int numcols, int parity)
1242 {
1243
1244         register jpc_fix_t *lptr;
1245         register jpc_fix_t *hptr;
1246         register int n;
1247         int llen;
1248
1249         llen = (numcols + 1 - parity) >> 1;
1250
1251         if (numcols > 1) {
1252
1253                 /* Apply the first lifting step. */
1254                 lptr = &a[0];
1255                 hptr = &a[llen];
1256                 if (!parity) {
1257                         lptr[0] -= (hptr[0] + 1) >> 1;
1258                         ++lptr;
1259                 }
1260                 n = llen - (!parity) - (parity != (numcols & 1));
1261                 while (n-- > 0) {
1262                         lptr[0] -= (hptr[0] + hptr[1] + 2) >> 2;
1263                         ++lptr;
1264                         ++hptr;
1265                 }
1266                 if (parity != (numcols & 1)) {
1267                         lptr[0] -= (hptr[0] + 1) >> 1;
1268                 }
1269
1270                 /* Apply the second lifting step. */
1271                 lptr = &a[0];
1272                 hptr = &a[llen];
1273                 if (parity) {
1274                         hptr[0] += lptr[0];
1275                         ++hptr;
1276                 }
1277                 n = numcols - llen - parity - (parity == (numcols & 1));
1278                 while (n-- > 0) {
1279                         hptr[0] += (lptr[0] + lptr[1]) >> 1;
1280                         ++hptr;
1281                         ++lptr;
1282                 }
1283                 if (parity == (numcols & 1)) {
1284                         hptr[0] += lptr[0];
1285                 }
1286
1287         } else {
1288
1289                 if (parity) {
1290                         lptr = &a[0];
1291                         lptr[0] >>= 1;
1292                 }
1293
1294         }
1295
1296 }
1297
1298 void jpc_ft_invlift_col(jpc_fix_t *a, int numrows, int stride, int parity)
1299 {
1300
1301         jpc_fix_t *lptr;
1302         jpc_fix_t *hptr;
1303 #if 0
1304         register jpc_fix_t *lptr2;
1305         register jpc_fix_t *hptr2;
1306         register int i;
1307 #endif
1308         register int n;
1309         int llen;
1310
1311         llen = (numrows + 1 - parity) >> 1;
1312
1313         if (numrows > 1) {
1314
1315                 /* Apply the first lifting step. */
1316                 lptr = &a[0];
1317                 hptr = &a[llen * stride];
1318                 if (!parity) {
1319                         lptr[0] -= (hptr[0] + 1) >> 1;
1320                         lptr += stride;
1321                 }
1322                 n = llen - (!parity) - (parity != (numrows & 1));
1323                 while (n-- > 0) {
1324                         lptr[0] -= (hptr[0] + hptr[stride] + 2) >> 2;
1325                         lptr += stride;
1326                         hptr += stride;
1327                 }
1328                 if (parity != (numrows & 1)) {
1329                         lptr[0] -= (hptr[0] + 1) >> 1;
1330                 }
1331
1332                 /* Apply the second lifting step. */
1333                 lptr = &a[0];
1334                 hptr = &a[llen * stride];
1335                 if (parity) {
1336                         hptr[0] += lptr[0];
1337                         hptr += stride;
1338                 }
1339                 n = numrows - llen - parity - (parity == (numrows & 1));
1340                 while (n-- > 0) {
1341                         hptr[0] += (lptr[0] + lptr[stride]) >> 1;
1342                         hptr += stride;
1343                         lptr += stride;
1344                 }
1345                 if (parity == (numrows & 1)) {
1346                         hptr[0] += lptr[0];
1347                 }
1348
1349         } else {
1350
1351                 if (parity) {
1352                         lptr = &a[0];
1353                         lptr[0] >>= 1;
1354                 }
1355
1356         }
1357
1358 }
1359
1360 void jpc_ft_invlift_colgrp(jpc_fix_t *a, int numrows, int stride, int parity)
1361 {
1362
1363         jpc_fix_t *lptr;
1364         jpc_fix_t *hptr;
1365         register jpc_fix_t *lptr2;
1366         register jpc_fix_t *hptr2;
1367         register int n;
1368         register int i;
1369         int llen;
1370
1371         llen = (numrows + 1 - parity) >> 1;
1372
1373         if (numrows > 1) {
1374
1375                 /* Apply the first lifting step. */
1376                 lptr = &a[0];
1377                 hptr = &a[llen * stride];
1378                 if (!parity) {
1379                         lptr2 = lptr;
1380                         hptr2 = hptr;
1381                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1382                                 lptr2[0] -= (hptr2[0] + 1) >> 1;
1383                                 ++lptr2;
1384                                 ++hptr2;
1385                         }
1386                         lptr += stride;
1387                 }
1388                 n = llen - (!parity) - (parity != (numrows & 1));
1389                 while (n-- > 0) {
1390                         lptr2 = lptr;
1391                         hptr2 = hptr;
1392                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1393                                 lptr2[0] -= (hptr2[0] + hptr2[stride] + 2) >> 2;
1394                                 ++lptr2;
1395                                 ++hptr2;
1396                         }
1397                         lptr += stride;
1398                         hptr += stride;
1399                 }
1400                 if (parity != (numrows & 1)) {
1401                         lptr2 = lptr;
1402                         hptr2 = hptr;
1403                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1404                                 lptr2[0] -= (hptr2[0] + 1) >> 1;
1405                                 ++lptr2;
1406                                 ++hptr2;
1407                         }
1408                 }
1409
1410                 /* Apply the second lifting step. */
1411                 lptr = &a[0];
1412                 hptr = &a[llen * stride];
1413                 if (parity) {
1414                         lptr2 = lptr;
1415                         hptr2 = hptr;
1416                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1417                                 hptr2[0] += lptr2[0];
1418                                 ++hptr2;
1419                                 ++lptr2;
1420                         }
1421                         hptr += stride;
1422                 }
1423                 n = numrows - llen - parity - (parity == (numrows & 1));
1424                 while (n-- > 0) {
1425                         lptr2 = lptr;
1426                         hptr2 = hptr;
1427                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1428                                 hptr2[0] += (lptr2[0] + lptr2[stride]) >> 1;
1429                                 ++lptr2;
1430                                 ++hptr2;
1431                         }
1432                         hptr += stride;
1433                         lptr += stride;
1434                 }
1435                 if (parity == (numrows & 1)) {
1436                         lptr2 = lptr;
1437                         hptr2 = hptr;
1438                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1439                                 hptr2[0] += lptr2[0];
1440                                 ++lptr2;
1441                                 ++hptr2;
1442                         }
1443                 }
1444
1445         } else {
1446
1447                 if (parity) {
1448                         lptr2 = &a[0];
1449                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1450                                 lptr2[0] >>= 1;
1451                                 ++lptr2;
1452                         }
1453                 }
1454
1455         }
1456
1457 }
1458
1459 void jpc_ft_invlift_colres(jpc_fix_t *a, int numrows, int numcols, int stride,
1460   int parity)
1461 {
1462
1463         jpc_fix_t *lptr;
1464         jpc_fix_t *hptr;
1465         register jpc_fix_t *lptr2;
1466         register jpc_fix_t *hptr2;
1467         register int n;
1468         register int i;
1469         int llen;
1470
1471         llen = (numrows + 1 - parity) >> 1;
1472
1473         if (numrows > 1) {
1474
1475                 /* Apply the first lifting step. */
1476                 lptr = &a[0];
1477                 hptr = &a[llen * stride];
1478                 if (!parity) {
1479                         lptr2 = lptr;
1480                         hptr2 = hptr;
1481                         for (i = 0; i < numcols; ++i) {
1482                                 lptr2[0] -= (hptr2[0] + 1) >> 1;
1483                                 ++lptr2;
1484                                 ++hptr2;
1485                         }
1486                         lptr += stride;
1487                 }
1488                 n = llen - (!parity) - (parity != (numrows & 1));
1489                 while (n-- > 0) {
1490                         lptr2 = lptr;
1491                         hptr2 = hptr;
1492                         for (i = 0; i < numcols; ++i) {
1493                                 lptr2[0] -= (hptr2[0] + hptr2[stride] + 2) >> 2;
1494                                 ++lptr2;
1495                                 ++hptr2;
1496                         }
1497                         lptr += stride;
1498                         hptr += stride;
1499                 }
1500                 if (parity != (numrows & 1)) {
1501                         lptr2 = lptr;
1502                         hptr2 = hptr;
1503                         for (i = 0; i < numcols; ++i) {
1504                                 lptr2[0] -= (hptr2[0] + 1) >> 1;
1505                                 ++lptr2;
1506                                 ++hptr2;
1507                         }
1508                 }
1509
1510                 /* Apply the second lifting step. */
1511                 lptr = &a[0];
1512                 hptr = &a[llen * stride];
1513                 if (parity) {
1514                         lptr2 = lptr;
1515                         hptr2 = hptr;
1516                         for (i = 0; i < numcols; ++i) {
1517                                 hptr2[0] += lptr2[0];
1518                                 ++hptr2;
1519                                 ++lptr2;
1520                         }
1521                         hptr += stride;
1522                 }
1523                 n = numrows - llen - parity - (parity == (numrows & 1));
1524                 while (n-- > 0) {
1525                         lptr2 = lptr;
1526                         hptr2 = hptr;
1527                         for (i = 0; i < numcols; ++i) {
1528                                 hptr2[0] += (lptr2[0] + lptr2[stride]) >> 1;
1529                                 ++lptr2;
1530                                 ++hptr2;
1531                         }
1532                         hptr += stride;
1533                         lptr += stride;
1534                 }
1535                 if (parity == (numrows & 1)) {
1536                         lptr2 = lptr;
1537                         hptr2 = hptr;
1538                         for (i = 0; i < numcols; ++i) {
1539                                 hptr2[0] += lptr2[0];
1540                                 ++lptr2;
1541                                 ++hptr2;
1542                         }
1543                 }
1544
1545         } else {
1546
1547                 if (parity) {
1548                         lptr2 = &a[0];
1549                         for (i = 0; i < numcols; ++i) {
1550                                 lptr2[0] >>= 1;
1551                                 ++lptr2;
1552                         }
1553                 }
1554
1555         }
1556
1557 }
1558
1559 int jpc_ft_analyze(jpc_fix_t *a, int xstart, int ystart, int width, int height,
1560   int stride)
1561 {
1562         int numrows = height;
1563         int numcols = width;
1564         int rowparity = ystart & 1;
1565         int colparity = xstart & 1;
1566         int i;
1567         jpc_fix_t *startptr;
1568         int maxcols;
1569
1570         maxcols = (numcols / JPC_QMFB_COLGRPSIZE) * JPC_QMFB_COLGRPSIZE;
1571         startptr = &a[0];
1572         for (i = 0; i < maxcols; i += JPC_QMFB_COLGRPSIZE) {
1573                 jpc_qmfb_split_colgrp(startptr, numrows, stride, rowparity);
1574                 jpc_ft_fwdlift_colgrp(startptr, numrows, stride, rowparity);
1575                 startptr += JPC_QMFB_COLGRPSIZE;
1576         }
1577         if (maxcols < numcols) {
1578                 jpc_qmfb_split_colres(startptr, numrows, numcols - maxcols, stride,
1579                   rowparity);
1580                 jpc_ft_fwdlift_colres(startptr, numrows, numcols - maxcols, stride,
1581                   rowparity);
1582         }
1583
1584         startptr = &a[0];
1585         for (i = 0; i < numrows; ++i) {
1586                 jpc_qmfb_split_row(startptr, numcols, colparity);
1587                 jpc_ft_fwdlift_row(startptr, numcols, colparity);
1588                 startptr += stride;
1589         }
1590
1591         return 0;
1592
1593 }
1594
1595 int jpc_ft_synthesize(int *a, int xstart, int ystart, int width, int height,
1596   int stride)
1597 {
1598         int numrows = height;
1599         int numcols = width;
1600         int rowparity = ystart & 1;
1601         int colparity = xstart & 1;
1602
1603         int maxcols;
1604         jpc_fix_t *startptr;
1605         int i;
1606
1607         startptr = &a[0];
1608         for (i = 0; i < numrows; ++i) {
1609                 jpc_ft_invlift_row(startptr, numcols, colparity);
1610                 jpc_qmfb_join_row(startptr, numcols, colparity);
1611                 startptr += stride;
1612         }
1613
1614         maxcols = (numcols / JPC_QMFB_COLGRPSIZE) * JPC_QMFB_COLGRPSIZE;
1615         startptr = &a[0];
1616         for (i = 0; i < maxcols; i += JPC_QMFB_COLGRPSIZE) {
1617                 jpc_ft_invlift_colgrp(startptr, numrows, stride, rowparity);
1618                 jpc_qmfb_join_colgrp(startptr, numrows, stride, rowparity);
1619                 startptr += JPC_QMFB_COLGRPSIZE;
1620         }
1621         if (maxcols < numcols) {
1622                 jpc_ft_invlift_colres(startptr, numrows, numcols - maxcols, stride,
1623                   rowparity);
1624                 jpc_qmfb_join_colres(startptr, numrows, numcols - maxcols, stride,
1625                   rowparity);
1626         }
1627
1628         return 0;
1629
1630 }
1631
1632 /******************************************************************************\
1633 * 9/7 transform
1634 \******************************************************************************/
1635
1636 #define ALPHA (-1.586134342059924)
1637 #define BETA (-0.052980118572961)
1638 #define GAMMA (0.882911075530934)
1639 #define DELTA (0.443506852043971)
1640 #define LGAIN (1.0 / 1.23017410558578)
1641 #define HGAIN (1.0 / 1.62578613134411)
1642
1643 void jpc_ns_fwdlift_row(jpc_fix_t *a, int numcols, int parity)
1644 {
1645
1646         register jpc_fix_t *lptr;
1647         register jpc_fix_t *hptr;
1648         register int n;
1649         int llen;
1650
1651         llen = (numcols + 1 - parity) >> 1;
1652
1653         if (numcols > 1) {
1654
1655                 /* Apply the first lifting step. */
1656                 lptr = &a[0];
1657                 hptr = &a[llen];
1658                 if (parity) {
1659                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
1660                           lptr[0]));
1661                         ++hptr;
1662                 }
1663                 n = numcols - llen - parity - (parity == (numcols & 1));
1664                 while (n-- > 0) {
1665                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
1666                           jpc_fix_add(lptr[0], lptr[1])));
1667                         ++hptr;
1668                         ++lptr;
1669                 }
1670                 if (parity == (numcols & 1)) {
1671                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
1672                           lptr[0]));
1673                 }
1674
1675                 /* Apply the second lifting step. */
1676                 lptr = &a[0];
1677                 hptr = &a[llen];
1678                 if (!parity) {
1679                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
1680                           hptr[0]));
1681                         ++lptr;
1682                 }
1683                 n = llen - (!parity) - (parity != (numcols & 1));
1684                 while (n-- > 0) {
1685                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(BETA),
1686                           jpc_fix_add(hptr[0], hptr[1])));
1687                         ++lptr;
1688                         ++hptr;
1689                 }
1690                 if (parity != (numcols & 1)) {
1691                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
1692                           hptr[0]));
1693                 }
1694
1695                 /* Apply the third lifting step. */
1696                 lptr = &a[0];
1697                 hptr = &a[llen];
1698                 if (parity) {
1699                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
1700                           lptr[0]));
1701                         ++hptr;
1702                 }
1703                 n = numcols - llen - parity - (parity == (numcols & 1));
1704                 while (n-- > 0) {
1705                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
1706                           jpc_fix_add(lptr[0], lptr[1])));
1707                         ++hptr;
1708                         ++lptr;
1709                 }
1710                 if (parity == (numcols & 1)) {
1711                         jpc_fix_pluseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
1712                           lptr[0]));
1713                 }
1714
1715                 /* Apply the fourth lifting step. */
1716                 lptr = &a[0];
1717                 hptr = &a[llen];
1718                 if (!parity) {
1719                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
1720                           hptr[0]));
1721                         ++lptr;
1722                 }
1723                 n = llen - (!parity) - (parity != (numcols & 1));
1724                 while (n-- > 0) {
1725                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(DELTA),
1726                           jpc_fix_add(hptr[0], hptr[1])));
1727                         ++lptr;
1728                         ++hptr;
1729                 }
1730                 if (parity != (numcols & 1)) {
1731                         jpc_fix_pluseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
1732                           hptr[0]));
1733                 }
1734
1735                 /* Apply the scaling step. */
1736 #if defined(WT_DOSCALE)
1737                 lptr = &a[0];
1738                 n = llen;
1739                 while (n-- > 0) {
1740                         lptr[0] = jpc_fix_mul(lptr[0], jpc_dbltofix(LGAIN));
1741                         ++lptr;
1742                 }
1743                 hptr = &a[llen];
1744                 n = numcols - llen;
1745                 while (n-- > 0) {
1746                         hptr[0] = jpc_fix_mul(hptr[0], jpc_dbltofix(HGAIN));
1747                         ++hptr;
1748                 }
1749 #endif
1750
1751         } else {
1752
1753 #if defined(WT_LENONE)
1754                 if (parity) {
1755                         lptr = &a[0];
1756                         lptr[0] <<= 1;
1757                 }
1758 #endif
1759
1760         }
1761
1762 }
1763
1764 void jpc_ns_fwdlift_colgrp(jpc_fix_t *a, int numrows, int stride,
1765   int parity)
1766 {
1767
1768         jpc_fix_t *lptr;
1769         jpc_fix_t *hptr;
1770         register jpc_fix_t *lptr2;
1771         register jpc_fix_t *hptr2;
1772         register int n;
1773         register int i;
1774         int llen;
1775
1776         llen = (numrows + 1 - parity) >> 1;
1777
1778         if (numrows > 1) {
1779
1780                 /* Apply the first lifting step. */
1781                 lptr = &a[0];
1782                 hptr = &a[llen * stride];
1783                 if (parity) {
1784                         lptr2 = lptr;
1785                         hptr2 = hptr;
1786                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1787                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
1788                                   lptr2[0]));
1789                                 ++hptr2;
1790                                 ++lptr2;
1791                         }
1792                         hptr += stride;
1793                 }
1794                 n = numrows - llen - parity - (parity == (numrows & 1));
1795                 while (n-- > 0) {
1796                         lptr2 = lptr;
1797                         hptr2 = hptr;
1798                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1799                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
1800                                   jpc_fix_add(lptr2[0], lptr2[stride])));
1801                                 ++lptr2;
1802                                 ++hptr2;
1803                         }
1804                         hptr += stride;
1805                         lptr += stride;
1806                 }
1807                 if (parity == (numrows & 1)) {
1808                         lptr2 = lptr;
1809                         hptr2 = hptr;
1810                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1811                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
1812                                   lptr2[0]));
1813                                 ++lptr2;
1814                                 ++hptr2;
1815                         }
1816                 }
1817
1818                 /* Apply the second lifting step. */
1819                 lptr = &a[0];
1820                 hptr = &a[llen * stride];
1821                 if (!parity) {
1822                         lptr2 = lptr;
1823                         hptr2 = hptr;
1824                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1825                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
1826                                   hptr2[0]));
1827                                 ++lptr2;
1828                                 ++hptr2;
1829                         }
1830                         lptr += stride;
1831                 }
1832                 n = llen - (!parity) - (parity != (numrows & 1));
1833                 while (n-- > 0) {
1834                         lptr2 = lptr;
1835                         hptr2 = hptr;
1836                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1837                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
1838                                   jpc_fix_add(hptr2[0], hptr2[stride])));
1839                                 ++lptr2;
1840                                 ++hptr2;
1841                         }
1842                         lptr += stride;
1843                         hptr += stride;
1844                 }
1845                 if (parity != (numrows & 1)) {
1846                         lptr2 = lptr;
1847                         hptr2 = hptr;
1848                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1849                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
1850                                   hptr2[0]));
1851                                 ++lptr2;
1852                                 ++hptr2;
1853                         }
1854                 }
1855
1856                 /* Apply the third lifting step. */
1857                 lptr = &a[0];
1858                 hptr = &a[llen * stride];
1859                 if (parity) {
1860                         lptr2 = lptr;
1861                         hptr2 = hptr;
1862                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1863                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
1864                                   lptr2[0]));
1865                                 ++hptr2;
1866                                 ++lptr2;
1867                         }
1868                         hptr += stride;
1869                 }
1870                 n = numrows - llen - parity - (parity == (numrows & 1));
1871                 while (n-- > 0) {
1872                         lptr2 = lptr;
1873                         hptr2 = hptr;
1874                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1875                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
1876                                   jpc_fix_add(lptr2[0], lptr2[stride])));
1877                                 ++lptr2;
1878                                 ++hptr2;
1879                         }
1880                         hptr += stride;
1881                         lptr += stride;
1882                 }
1883                 if (parity == (numrows & 1)) {
1884                         lptr2 = lptr;
1885                         hptr2 = hptr;
1886                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1887                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
1888                                   lptr2[0]));
1889                                 ++lptr2;
1890                                 ++hptr2;
1891                         }
1892                 }
1893
1894                 /* Apply the fourth lifting step. */
1895                 lptr = &a[0];
1896                 hptr = &a[llen * stride];
1897                 if (!parity) {
1898                         lptr2 = lptr;
1899                         hptr2 = hptr;
1900                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1901                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
1902                                   hptr2[0]));
1903                                 ++lptr2;
1904                                 ++hptr2;
1905                         }
1906                         lptr += stride;
1907                 }
1908                 n = llen - (!parity) - (parity != (numrows & 1));
1909                 while (n-- > 0) {
1910                         lptr2 = lptr;
1911                         hptr2 = hptr;
1912                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1913                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
1914                                   jpc_fix_add(hptr2[0], hptr2[stride])));
1915                                 ++lptr2;
1916                                 ++hptr2;
1917                         }
1918                         lptr += stride;
1919                         hptr += stride;
1920                 }
1921                 if (parity != (numrows & 1)) {
1922                         lptr2 = lptr;
1923                         hptr2 = hptr;
1924                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1925                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
1926                                   hptr2[0]));
1927                                 ++lptr2;
1928                                 ++hptr2;
1929                         }
1930                 }
1931
1932                 /* Apply the scaling step. */
1933 #if defined(WT_DOSCALE)
1934                 lptr = &a[0];
1935                 n = llen;
1936                 while (n-- > 0) {
1937                         lptr2 = lptr;
1938                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1939                                 lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(LGAIN));
1940                                 ++lptr2;
1941                         }
1942                         lptr += stride;
1943                 }
1944                 hptr = &a[llen * stride];
1945                 n = numrows - llen;
1946                 while (n-- > 0) {
1947                         hptr2 = hptr;
1948                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1949                                 hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(HGAIN));
1950                                 ++hptr2;
1951                         }
1952                         hptr += stride;
1953                 }
1954 #endif
1955
1956         } else {
1957
1958 #if defined(WT_LENONE)
1959                 if (parity) {
1960                         lptr2 = &a[0];
1961                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
1962                                 lptr2[0] <<= 1;
1963                                 ++lptr2;
1964                         }
1965                 }
1966 #endif
1967
1968         }
1969
1970 }
1971
1972 void jpc_ns_fwdlift_colres(jpc_fix_t *a, int numrows, int numcols,
1973   int stride, int parity)
1974 {
1975
1976         jpc_fix_t *lptr;
1977         jpc_fix_t *hptr;
1978         register jpc_fix_t *lptr2;
1979         register jpc_fix_t *hptr2;
1980         register int n;
1981         register int i;
1982         int llen;
1983
1984         llen = (numrows + 1 - parity) >> 1;
1985
1986         if (numrows > 1) {
1987
1988                 /* Apply the first lifting step. */
1989                 lptr = &a[0];
1990                 hptr = &a[llen * stride];
1991                 if (parity) {
1992                         lptr2 = lptr;
1993                         hptr2 = hptr;
1994                         for (i = 0; i < numcols; ++i) {
1995                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
1996                                   lptr2[0]));
1997                                 ++hptr2;
1998                                 ++lptr2;
1999                         }
2000                         hptr += stride;
2001                 }
2002                 n = numrows - llen - parity - (parity == (numrows & 1));
2003                 while (n-- > 0) {
2004                         lptr2 = lptr;
2005                         hptr2 = hptr;
2006                         for (i = 0; i < numcols; ++i) {
2007                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
2008                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2009                                 ++lptr2;
2010                                 ++hptr2;
2011                         }
2012                         hptr += stride;
2013                         lptr += stride;
2014                 }
2015                 if (parity == (numrows & 1)) {
2016                         lptr2 = lptr;
2017                         hptr2 = hptr;
2018                         for (i = 0; i < numcols; ++i) {
2019                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
2020                                   lptr2[0]));
2021                                 ++lptr2;
2022                                 ++hptr2;
2023                         }
2024                 }
2025
2026                 /* Apply the second lifting step. */
2027                 lptr = &a[0];
2028                 hptr = &a[llen * stride];
2029                 if (!parity) {
2030                         lptr2 = lptr;
2031                         hptr2 = hptr;
2032                         for (i = 0; i < numcols; ++i) {
2033                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2034                                   hptr2[0]));
2035                                 ++lptr2;
2036                                 ++hptr2;
2037                         }
2038                         lptr += stride;
2039                 }
2040                 n = llen - (!parity) - (parity != (numrows & 1));
2041                 while (n-- > 0) {
2042                         lptr2 = lptr;
2043                         hptr2 = hptr;
2044                         for (i = 0; i < numcols; ++i) {
2045                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
2046                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2047                                 ++lptr2;
2048                                 ++hptr2;
2049                         }
2050                         lptr += stride;
2051                         hptr += stride;
2052                 }
2053                 if (parity != (numrows & 1)) {
2054                         lptr2 = lptr;
2055                         hptr2 = hptr;
2056                         for (i = 0; i < numcols; ++i) {
2057                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2058                                   hptr2[0]));
2059                                 ++lptr2;
2060                                 ++hptr2;
2061                         }
2062                 }
2063
2064                 /* Apply the third lifting step. */
2065                 lptr = &a[0];
2066                 hptr = &a[llen * stride];
2067                 if (parity) {
2068                         lptr2 = lptr;
2069                         hptr2 = hptr;
2070                         for (i = 0; i < numcols; ++i) {
2071                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2072                                   lptr2[0]));
2073                                 ++hptr2;
2074                                 ++lptr2;
2075                         }
2076                         hptr += stride;
2077                 }
2078                 n = numrows - llen - parity - (parity == (numrows & 1));
2079                 while (n-- > 0) {
2080                         lptr2 = lptr;
2081                         hptr2 = hptr;
2082                         for (i = 0; i < numcols; ++i) {
2083                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2084                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2085                                 ++lptr2;
2086                                 ++hptr2;
2087                         }
2088                         hptr += stride;
2089                         lptr += stride;
2090                 }
2091                 if (parity == (numrows & 1)) {
2092                         lptr2 = lptr;
2093                         hptr2 = hptr;
2094                         for (i = 0; i < numcols; ++i) {
2095                                 jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2096                                   lptr2[0]));
2097                                 ++lptr2;
2098                                 ++hptr2;
2099                         }
2100                 }
2101
2102                 /* Apply the fourth lifting step. */
2103                 lptr = &a[0];
2104                 hptr = &a[llen * stride];
2105                 if (!parity) {
2106                         lptr2 = lptr;
2107                         hptr2 = hptr;
2108                         for (i = 0; i < numcols; ++i) {
2109                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2110                                   hptr2[0]));
2111                                 ++lptr2;
2112                                 ++hptr2;
2113                         }
2114                         lptr += stride;
2115                 }
2116                 n = llen - (!parity) - (parity != (numrows & 1));
2117                 while (n-- > 0) {
2118                         lptr2 = lptr;
2119                         hptr2 = hptr;
2120                         for (i = 0; i < numcols; ++i) {
2121                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2122                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2123                                 ++lptr2;
2124                                 ++hptr2;
2125                         }
2126                         lptr += stride;
2127                         hptr += stride;
2128                 }
2129                 if (parity != (numrows & 1)) {
2130                         lptr2 = lptr;
2131                         hptr2 = hptr;
2132                         for (i = 0; i < numcols; ++i) {
2133                                 jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2134                                   hptr2[0]));
2135                                 ++lptr2;
2136                                 ++hptr2;
2137                         }
2138                 }
2139
2140                 /* Apply the scaling step. */
2141 #if defined(WT_DOSCALE)
2142                 lptr = &a[0];
2143                 n = llen;
2144                 while (n-- > 0) {
2145                         lptr2 = lptr;
2146                         for (i = 0; i < numcols; ++i) {
2147                                 lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(LGAIN));
2148                                 ++lptr2;
2149                         }
2150                         lptr += stride;
2151                 }
2152                 hptr = &a[llen * stride];
2153                 n = numrows - llen;
2154                 while (n-- > 0) {
2155                         hptr2 = hptr;
2156                         for (i = 0; i < numcols; ++i) {
2157                                 hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(HGAIN));
2158                                 ++hptr2;
2159                         }
2160                         hptr += stride;
2161                 }
2162 #endif
2163
2164         } else {
2165
2166 #if defined(WT_LENONE)
2167                 if (parity) {
2168                         lptr2 = &a[0];
2169                         for (i = 0; i < numcols; ++i) {
2170                                 lptr2[0] <<= 1;
2171                                 ++lptr2;
2172                         }
2173                 }
2174 #endif
2175
2176         }
2177
2178 }
2179
2180 void jpc_ns_fwdlift_col(jpc_fix_t *a, int numrows, int stride,
2181   int parity)
2182 {
2183
2184         jpc_fix_t *lptr;
2185         jpc_fix_t *hptr;
2186         register jpc_fix_t *lptr2;
2187         register jpc_fix_t *hptr2;
2188         register int n;
2189         int llen;
2190
2191         llen = (numrows + 1 - parity) >> 1;
2192
2193         if (numrows > 1) {
2194
2195                 /* Apply the first lifting step. */
2196                 lptr = &a[0];
2197                 hptr = &a[llen * stride];
2198                 if (parity) {
2199                         lptr2 = lptr;
2200                         hptr2 = hptr;
2201                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
2202                           lptr2[0]));
2203                         ++hptr2;
2204                         ++lptr2;
2205                         hptr += stride;
2206                 }
2207                 n = numrows - llen - parity - (parity == (numrows & 1));
2208                 while (n-- > 0) {
2209                         lptr2 = lptr;
2210                         hptr2 = hptr;
2211                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
2212                           jpc_fix_add(lptr2[0], lptr2[stride])));
2213                         ++lptr2;
2214                         ++hptr2;
2215                         hptr += stride;
2216                         lptr += stride;
2217                 }
2218                 if (parity == (numrows & 1)) {
2219                         lptr2 = lptr;
2220                         hptr2 = hptr;
2221                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
2222                           lptr2[0]));
2223                         ++lptr2;
2224                         ++hptr2;
2225                 }
2226
2227                 /* Apply the second lifting step. */
2228                 lptr = &a[0];
2229                 hptr = &a[llen * stride];
2230                 if (!parity) {
2231                         lptr2 = lptr;
2232                         hptr2 = hptr;
2233                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2234                           hptr2[0]));
2235                         ++lptr2;
2236                         ++hptr2;
2237                         lptr += stride;
2238                 }
2239                 n = llen - (!parity) - (parity != (numrows & 1));
2240                 while (n-- > 0) {
2241                         lptr2 = lptr;
2242                         hptr2 = hptr;
2243                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
2244                           jpc_fix_add(hptr2[0], hptr2[stride])));
2245                         ++lptr2;
2246                         ++hptr2;
2247                         lptr += stride;
2248                         hptr += stride;
2249                 }
2250                 if (parity != (numrows & 1)) {
2251                         lptr2 = lptr;
2252                         hptr2 = hptr;
2253                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2254                           hptr2[0]));
2255                         ++lptr2;
2256                         ++hptr2;
2257                 }
2258
2259                 /* Apply the third lifting step. */
2260                 lptr = &a[0];
2261                 hptr = &a[llen * stride];
2262                 if (parity) {
2263                         lptr2 = lptr;
2264                         hptr2 = hptr;
2265                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2266                           lptr2[0]));
2267                         ++hptr2;
2268                         ++lptr2;
2269                         hptr += stride;
2270                 }
2271                 n = numrows - llen - parity - (parity == (numrows & 1));
2272                 while (n-- > 0) {
2273                         lptr2 = lptr;
2274                         hptr2 = hptr;
2275                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2276                           jpc_fix_add(lptr2[0], lptr2[stride])));
2277                         ++lptr2;
2278                         ++hptr2;
2279                         hptr += stride;
2280                         lptr += stride;
2281                 }
2282                 if (parity == (numrows & 1)) {
2283                         lptr2 = lptr;
2284                         hptr2 = hptr;
2285                         jpc_fix_pluseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2286                           lptr2[0]));
2287                         ++lptr2;
2288                         ++hptr2;
2289                 }
2290
2291                 /* Apply the fourth lifting step. */
2292                 lptr = &a[0];
2293                 hptr = &a[llen * stride];
2294                 if (!parity) {
2295                         lptr2 = lptr;
2296                         hptr2 = hptr;
2297                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2298                           hptr2[0]));
2299                         ++lptr2;
2300                         ++hptr2;
2301                         lptr += stride;
2302                 }
2303                 n = llen - (!parity) - (parity != (numrows & 1));
2304                 while (n-- > 0) {
2305                         lptr2 = lptr;
2306                         hptr2 = hptr;
2307                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2308                           jpc_fix_add(hptr2[0], hptr2[stride])));
2309                         ++lptr2;
2310                         ++hptr2;
2311                         lptr += stride;
2312                         hptr += stride;
2313                 }
2314                 if (parity != (numrows & 1)) {
2315                         lptr2 = lptr;
2316                         hptr2 = hptr;
2317                         jpc_fix_pluseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2318                           hptr2[0]));
2319                         ++lptr2;
2320                         ++hptr2;
2321                 }
2322
2323                 /* Apply the scaling step. */
2324 #if defined(WT_DOSCALE)
2325                 lptr = &a[0];
2326                 n = llen;
2327                 while (n-- > 0) {
2328                         lptr2 = lptr;
2329                         lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(LGAIN));
2330                         ++lptr2;
2331                         lptr += stride;
2332                 }
2333                 hptr = &a[llen * stride];
2334                 n = numrows - llen;
2335                 while (n-- > 0) {
2336                         hptr2 = hptr;
2337                         hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(HGAIN));
2338                         ++hptr2;
2339                         hptr += stride;
2340                 }
2341 #endif
2342
2343         } else {
2344
2345 #if defined(WT_LENONE)
2346                 if (parity) {
2347                         lptr2 = &a[0];
2348                         lptr2[0] <<= 1;
2349                         ++lptr2;
2350                 }
2351 #endif
2352
2353         }
2354
2355 }
2356
2357 void jpc_ns_invlift_row(jpc_fix_t *a, int numcols, int parity)
2358 {
2359
2360         register jpc_fix_t *lptr;
2361         register jpc_fix_t *hptr;
2362         register int n;
2363         int llen;
2364
2365         llen = (numcols + 1 - parity) >> 1;
2366
2367         if (numcols > 1) {
2368
2369                 /* Apply the scaling step. */
2370 #if defined(WT_DOSCALE)
2371                 lptr = &a[0];
2372                 n = llen;
2373                 while (n-- > 0) {
2374                         lptr[0] = jpc_fix_mul(lptr[0], jpc_dbltofix(1.0 / LGAIN));
2375                         ++lptr;
2376                 }
2377                 hptr = &a[llen];
2378                 n = numcols - llen;
2379                 while (n-- > 0) {
2380                         hptr[0] = jpc_fix_mul(hptr[0], jpc_dbltofix(1.0 / HGAIN));
2381                         ++hptr;
2382                 }
2383 #endif
2384
2385                 /* Apply the first lifting step. */
2386                 lptr = &a[0];
2387                 hptr = &a[llen];
2388                 if (!parity) {
2389                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2390                           hptr[0]));
2391                         ++lptr;
2392                 }
2393                 n = llen - (!parity) - (parity != (numcols & 1));
2394                 while (n-- > 0) {
2395                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2396                           jpc_fix_add(hptr[0], hptr[1])));
2397                         ++lptr;
2398                         ++hptr;
2399                 }
2400                 if (parity != (numcols & 1)) {
2401                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * DELTA),
2402                           hptr[0]));
2403                 }
2404
2405                 /* Apply the second lifting step. */
2406                 lptr = &a[0];
2407                 hptr = &a[llen];
2408                 if (parity) {
2409                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2410                           lptr[0]));
2411                         ++hptr;
2412                 }
2413                 n = numcols - llen - parity - (parity == (numcols & 1));
2414                 while (n-- > 0) {
2415                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2416                           jpc_fix_add(lptr[0], lptr[1])));
2417                         ++hptr;
2418                         ++lptr;
2419                 }
2420                 if (parity == (numcols & 1)) {
2421                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * GAMMA),
2422                           lptr[0]));
2423                 }
2424
2425                 /* Apply the third lifting step. */
2426                 lptr = &a[0];
2427                 hptr = &a[llen];
2428                 if (!parity) {
2429                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2430                           hptr[0]));
2431                         ++lptr;
2432                 }
2433                 n = llen - (!parity) - (parity != (numcols & 1));
2434                 while (n-- > 0) {
2435                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(BETA),
2436                           jpc_fix_add(hptr[0], hptr[1])));
2437                         ++lptr;
2438                         ++hptr;
2439                 }
2440                 if (parity != (numcols & 1)) {
2441                         jpc_fix_minuseq(lptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2442                           hptr[0]));
2443                 }
2444
2445                 /* Apply the fourth lifting step. */
2446                 lptr = &a[0];
2447                 hptr = &a[llen];
2448                 if (parity) {
2449                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
2450                           lptr[0]));
2451                         ++hptr;
2452                 }
2453                 n = numcols - llen - parity - (parity == (numcols & 1));
2454                 while (n-- > 0) {
2455                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
2456                           jpc_fix_add(lptr[0], lptr[1])));
2457                         ++hptr;
2458                         ++lptr;
2459                 }
2460                 if (parity == (numcols & 1)) {
2461                         jpc_fix_minuseq(hptr[0], jpc_fix_mul(jpc_dbltofix(2.0 * ALPHA),
2462                           lptr[0]));
2463                 }
2464
2465         } else {
2466
2467 #if defined(WT_LENONE)
2468                 if (parity) {
2469                         lptr = &a[0];
2470                         lptr[0] >>= 1;
2471                 }
2472 #endif
2473
2474         }
2475
2476 }
2477
2478 void jpc_ns_invlift_colgrp(jpc_fix_t *a, int numrows, int stride,
2479   int parity)
2480 {
2481
2482         jpc_fix_t *lptr;
2483         jpc_fix_t *hptr;
2484         register jpc_fix_t *lptr2;
2485         register jpc_fix_t *hptr2;
2486         register int n;
2487         register int i;
2488         int llen;
2489
2490         llen = (numrows + 1 - parity) >> 1;
2491
2492         if (numrows > 1) {
2493
2494                 /* Apply the scaling step. */
2495 #if defined(WT_DOSCALE)
2496                 lptr = &a[0];
2497                 n = llen;
2498                 while (n-- > 0) {
2499                         lptr2 = lptr;
2500                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2501                                 lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(1.0 / LGAIN));
2502                                 ++lptr2;
2503                         }
2504                         lptr += stride;
2505                 }
2506                 hptr = &a[llen * stride];
2507                 n = numrows - llen;
2508                 while (n-- > 0) {
2509                         hptr2 = hptr;
2510                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2511                                 hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(1.0 / HGAIN));
2512                                 ++hptr2;
2513                         }
2514                         hptr += stride;
2515                 }
2516 #endif
2517
2518                 /* Apply the first lifting step. */
2519                 lptr = &a[0];
2520                 hptr = &a[llen * stride];
2521                 if (!parity) {
2522                         lptr2 = lptr;
2523                         hptr2 = hptr;
2524                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2525                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2526                                   DELTA), hptr2[0]));
2527                                 ++lptr2;
2528                                 ++hptr2;
2529                         }
2530                         lptr += stride;
2531                 }
2532                 n = llen - (!parity) - (parity != (numrows & 1));
2533                 while (n-- > 0) {
2534                         lptr2 = lptr;
2535                         hptr2 = hptr;
2536                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2537                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2538                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2539                                 ++lptr2;
2540                                 ++hptr2;
2541                         }
2542                         lptr += stride;
2543                         hptr += stride;
2544                 }
2545                 if (parity != (numrows & 1)) {
2546                         lptr2 = lptr;
2547                         hptr2 = hptr;
2548                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2549                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2550                                   DELTA), hptr2[0]));
2551                                 ++lptr2;
2552                                 ++hptr2;
2553                         }
2554                 }
2555
2556                 /* Apply the second lifting step. */
2557                 lptr = &a[0];
2558                 hptr = &a[llen * stride];
2559                 if (parity) {
2560                         lptr2 = lptr;
2561                         hptr2 = hptr;
2562                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2563                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2564                                   GAMMA), lptr2[0]));
2565                                 ++hptr2;
2566                                 ++lptr2;
2567                         }
2568                         hptr += stride;
2569                 }
2570                 n = numrows - llen - parity - (parity == (numrows & 1));
2571                 while (n-- > 0) {
2572                         lptr2 = lptr;
2573                         hptr2 = hptr;
2574                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2575                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2576                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2577                                 ++lptr2;
2578                                 ++hptr2;
2579                         }
2580                         hptr += stride;
2581                         lptr += stride;
2582                 }
2583                 if (parity == (numrows & 1)) {
2584                         lptr2 = lptr;
2585                         hptr2 = hptr;
2586                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2587                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2588                                   GAMMA), lptr2[0]));
2589                                 ++lptr2;
2590                                 ++hptr2;
2591                         }
2592                 }
2593
2594                 /* Apply the third lifting step. */
2595                 lptr = &a[0];
2596                 hptr = &a[llen * stride];
2597                 if (!parity) {
2598                         lptr2 = lptr;
2599                         hptr2 = hptr;
2600                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2601                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2602                                   hptr2[0]));
2603                                 ++lptr2;
2604                                 ++hptr2;
2605                         }
2606                         lptr += stride;
2607                 }
2608                 n = llen - (!parity) - (parity != (numrows & 1));
2609                 while (n-- > 0) {
2610                         lptr2 = lptr;
2611                         hptr2 = hptr;
2612                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2613                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
2614                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2615                                 ++lptr2;
2616                                 ++hptr2;
2617                         }
2618                         lptr += stride;
2619                         hptr += stride;
2620                 }
2621                 if (parity != (numrows & 1)) {
2622                         lptr2 = lptr;
2623                         hptr2 = hptr;
2624                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2625                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2626                                   hptr2[0]));
2627                                 ++lptr2;
2628                                 ++hptr2;
2629                         }
2630                 }
2631
2632                 /* Apply the fourth lifting step. */
2633                 lptr = &a[0];
2634                 hptr = &a[llen * stride];
2635                 if (parity) {
2636                         lptr2 = lptr;
2637                         hptr2 = hptr;
2638                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2639                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2640                                   ALPHA), lptr2[0]));
2641                                 ++hptr2;
2642                                 ++lptr2;
2643                         }
2644                         hptr += stride;
2645                 }
2646                 n = numrows - llen - parity - (parity == (numrows & 1));
2647                 while (n-- > 0) {
2648                         lptr2 = lptr;
2649                         hptr2 = hptr;
2650                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2651                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
2652                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2653                                 ++lptr2;
2654                                 ++hptr2;
2655                         }
2656                         hptr += stride;
2657                         lptr += stride;
2658                 }
2659                 if (parity == (numrows & 1)) {
2660                         lptr2 = lptr;
2661                         hptr2 = hptr;
2662                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2663                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2664                                   ALPHA), lptr2[0]));
2665                                 ++lptr2;
2666                                 ++hptr2;
2667                         }
2668                 }
2669
2670         } else {
2671
2672 #if defined(WT_LENONE)
2673                 if (parity) {
2674                         lptr2 = &a[0];
2675                         for (i = 0; i < JPC_QMFB_COLGRPSIZE; ++i) {
2676                                 lptr2[0] >>= 1;
2677                                 ++lptr2;
2678                         }
2679                 }
2680 #endif
2681
2682         }
2683
2684 }
2685
2686 void jpc_ns_invlift_colres(jpc_fix_t *a, int numrows, int numcols,
2687   int stride, int parity)
2688 {
2689
2690         jpc_fix_t *lptr;
2691         jpc_fix_t *hptr;
2692         register jpc_fix_t *lptr2;
2693         register jpc_fix_t *hptr2;
2694         register int n;
2695         register int i;
2696         int llen;
2697
2698         llen = (numrows + 1 - parity) >> 1;
2699
2700         if (numrows > 1) {
2701
2702                 /* Apply the scaling step. */
2703 #if defined(WT_DOSCALE)
2704                 lptr = &a[0];
2705                 n = llen;
2706                 while (n-- > 0) {
2707                         lptr2 = lptr;
2708                         for (i = 0; i < numcols; ++i) {
2709                                 lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(1.0 / LGAIN));
2710                                 ++lptr2;
2711                         }
2712                         lptr += stride;
2713                 }
2714                 hptr = &a[llen * stride];
2715                 n = numrows - llen;
2716                 while (n-- > 0) {
2717                         hptr2 = hptr;
2718                         for (i = 0; i < numcols; ++i) {
2719                                 hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(1.0 / HGAIN));
2720                                 ++hptr2;
2721                         }
2722                         hptr += stride;
2723                 }
2724 #endif
2725
2726                 /* Apply the first lifting step. */
2727                 lptr = &a[0];
2728                 hptr = &a[llen * stride];
2729                 if (!parity) {
2730                         lptr2 = lptr;
2731                         hptr2 = hptr;
2732                         for (i = 0; i < numcols; ++i) {
2733                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2734                                   DELTA), hptr2[0]));
2735                                 ++lptr2;
2736                                 ++hptr2;
2737                         }
2738                         lptr += stride;
2739                 }
2740                 n = llen - (!parity) - (parity != (numrows & 1));
2741                 while (n-- > 0) {
2742                         lptr2 = lptr;
2743                         hptr2 = hptr;
2744                         for (i = 0; i < numcols; ++i) {
2745                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2746                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2747                                 ++lptr2;
2748                                 ++hptr2;
2749                         }
2750                         lptr += stride;
2751                         hptr += stride;
2752                 }
2753                 if (parity != (numrows & 1)) {
2754                         lptr2 = lptr;
2755                         hptr2 = hptr;
2756                         for (i = 0; i < numcols; ++i) {
2757                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2758                                   DELTA), hptr2[0]));
2759                                 ++lptr2;
2760                                 ++hptr2;
2761                         }
2762                 }
2763
2764                 /* Apply the second lifting step. */
2765                 lptr = &a[0];
2766                 hptr = &a[llen * stride];
2767                 if (parity) {
2768                         lptr2 = lptr;
2769                         hptr2 = hptr;
2770                         for (i = 0; i < numcols; ++i) {
2771                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2772                                   GAMMA), lptr2[0]));
2773                                 ++hptr2;
2774                                 ++lptr2;
2775                         }
2776                         hptr += stride;
2777                 }
2778                 n = numrows - llen - parity - (parity == (numrows & 1));
2779                 while (n-- > 0) {
2780                         lptr2 = lptr;
2781                         hptr2 = hptr;
2782                         for (i = 0; i < numcols; ++i) {
2783                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2784                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2785                                 ++lptr2;
2786                                 ++hptr2;
2787                         }
2788                         hptr += stride;
2789                         lptr += stride;
2790                 }
2791                 if (parity == (numrows & 1)) {
2792                         lptr2 = lptr;
2793                         hptr2 = hptr;
2794                         for (i = 0; i < numcols; ++i) {
2795                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2796                                   GAMMA), lptr2[0]));
2797                                 ++lptr2;
2798                                 ++hptr2;
2799                         }
2800                 }
2801
2802                 /* Apply the third lifting step. */
2803                 lptr = &a[0];
2804                 hptr = &a[llen * stride];
2805                 if (!parity) {
2806                         lptr2 = lptr;
2807                         hptr2 = hptr;
2808                         for (i = 0; i < numcols; ++i) {
2809                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2810                                   hptr2[0]));
2811                                 ++lptr2;
2812                                 ++hptr2;
2813                         }
2814                         lptr += stride;
2815                 }
2816                 n = llen - (!parity) - (parity != (numrows & 1));
2817                 while (n-- > 0) {
2818                         lptr2 = lptr;
2819                         hptr2 = hptr;
2820                         for (i = 0; i < numcols; ++i) {
2821                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
2822                                   jpc_fix_add(hptr2[0], hptr2[stride])));
2823                                 ++lptr2;
2824                                 ++hptr2;
2825                         }
2826                         lptr += stride;
2827                         hptr += stride;
2828                 }
2829                 if (parity != (numrows & 1)) {
2830                         lptr2 = lptr;
2831                         hptr2 = hptr;
2832                         for (i = 0; i < numcols; ++i) {
2833                                 jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
2834                                   hptr2[0]));
2835                                 ++lptr2;
2836                                 ++hptr2;
2837                         }
2838                 }
2839
2840                 /* Apply the fourth lifting step. */
2841                 lptr = &a[0];
2842                 hptr = &a[llen * stride];
2843                 if (parity) {
2844                         lptr2 = lptr;
2845                         hptr2 = hptr;
2846                         for (i = 0; i < numcols; ++i) {
2847                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2848                                   ALPHA), lptr2[0]));
2849                                 ++hptr2;
2850                                 ++lptr2;
2851                         }
2852                         hptr += stride;
2853                 }
2854                 n = numrows - llen - parity - (parity == (numrows & 1));
2855                 while (n-- > 0) {
2856                         lptr2 = lptr;
2857                         hptr2 = hptr;
2858                         for (i = 0; i < numcols; ++i) {
2859                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
2860                                   jpc_fix_add(lptr2[0], lptr2[stride])));
2861                                 ++lptr2;
2862                                 ++hptr2;
2863                         }
2864                         hptr += stride;
2865                         lptr += stride;
2866                 }
2867                 if (parity == (numrows & 1)) {
2868                         lptr2 = lptr;
2869                         hptr2 = hptr;
2870                         for (i = 0; i < numcols; ++i) {
2871                                 jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2872                                   ALPHA), lptr2[0]));
2873                                 ++lptr2;
2874                                 ++hptr2;
2875                         }
2876                 }
2877
2878         } else {
2879
2880 #if defined(WT_LENONE)
2881                 if (parity) {
2882                         lptr2 = &a[0];
2883                         for (i = 0; i < numcols; ++i) {
2884                                 lptr2[0] >>= 1;
2885                                 ++lptr2;
2886                         }
2887                 }
2888 #endif
2889
2890         }
2891
2892 }
2893
2894 void jpc_ns_invlift_col(jpc_fix_t *a, int numrows, int stride,
2895   int parity)
2896 {
2897
2898         jpc_fix_t *lptr;
2899         jpc_fix_t *hptr;
2900         register jpc_fix_t *lptr2;
2901         register jpc_fix_t *hptr2;
2902         register int n;
2903         int llen;
2904
2905         llen = (numrows + 1 - parity) >> 1;
2906
2907         if (numrows > 1) {
2908
2909                 /* Apply the scaling step. */
2910 #if defined(WT_DOSCALE)
2911                 lptr = &a[0];
2912                 n = llen;
2913                 while (n-- > 0) {
2914                         lptr2 = lptr;
2915                         lptr2[0] = jpc_fix_mul(lptr2[0], jpc_dbltofix(1.0 / LGAIN));
2916                         ++lptr2;
2917                         lptr += stride;
2918                 }
2919                 hptr = &a[llen * stride];
2920                 n = numrows - llen;
2921                 while (n-- > 0) {
2922                         hptr2 = hptr;
2923                         hptr2[0] = jpc_fix_mul(hptr2[0], jpc_dbltofix(1.0 / HGAIN));
2924                         ++hptr2;
2925                         hptr += stride;
2926                 }
2927 #endif
2928
2929                 /* Apply the first lifting step. */
2930                 lptr = &a[0];
2931                 hptr = &a[llen * stride];
2932                 if (!parity) {
2933                         lptr2 = lptr;
2934                         hptr2 = hptr;
2935                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2936                           DELTA), hptr2[0]));
2937                         ++lptr2;
2938                         ++hptr2;
2939                         lptr += stride;
2940                 }
2941                 n = llen - (!parity) - (parity != (numrows & 1));
2942                 while (n-- > 0) {
2943                         lptr2 = lptr;
2944                         hptr2 = hptr;
2945                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(DELTA),
2946                           jpc_fix_add(hptr2[0], hptr2[stride])));
2947                         ++lptr2;
2948                         ++hptr2;
2949                         lptr += stride;
2950                         hptr += stride;
2951                 }
2952                 if (parity != (numrows & 1)) {
2953                         lptr2 = lptr;
2954                         hptr2 = hptr;
2955                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2956                           DELTA), hptr2[0]));
2957                         ++lptr2;
2958                         ++hptr2;
2959                 }
2960
2961                 /* Apply the second lifting step. */
2962                 lptr = &a[0];
2963                 hptr = &a[llen * stride];
2964                 if (parity) {
2965                         lptr2 = lptr;
2966                         hptr2 = hptr;
2967                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2968                           GAMMA), lptr2[0]));
2969                         ++hptr2;
2970                         ++lptr2;
2971                         hptr += stride;
2972                 }
2973                 n = numrows - llen - parity - (parity == (numrows & 1));
2974                 while (n-- > 0) {
2975                         lptr2 = lptr;
2976                         hptr2 = hptr;
2977                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(GAMMA),
2978                           jpc_fix_add(lptr2[0], lptr2[stride])));
2979                         ++lptr2;
2980                         ++hptr2;
2981                         hptr += stride;
2982                         lptr += stride;
2983                 }
2984                 if (parity == (numrows & 1)) {
2985                         lptr2 = lptr;
2986                         hptr2 = hptr;
2987                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
2988                           GAMMA), lptr2[0]));
2989                         ++lptr2;
2990                         ++hptr2;
2991                 }
2992
2993                 /* Apply the third lifting step. */
2994                 lptr = &a[0];
2995                 hptr = &a[llen * stride];
2996                 if (!parity) {
2997                         lptr2 = lptr;
2998                         hptr2 = hptr;
2999                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
3000                           hptr2[0]));
3001                         ++lptr2;
3002                         ++hptr2;
3003                         lptr += stride;
3004                 }
3005                 n = llen - (!parity) - (parity != (numrows & 1));
3006                 while (n-- > 0) {
3007                         lptr2 = lptr;
3008                         hptr2 = hptr;
3009                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(BETA),
3010                           jpc_fix_add(hptr2[0], hptr2[stride])));
3011                         ++lptr2;
3012                         ++hptr2;
3013                         lptr += stride;
3014                         hptr += stride;
3015                 }
3016                 if (parity != (numrows & 1)) {
3017                         lptr2 = lptr;
3018                         hptr2 = hptr;
3019                         jpc_fix_minuseq(lptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 * BETA),
3020                           hptr2[0]));
3021                         ++lptr2;
3022                         ++hptr2;
3023                 }
3024
3025                 /* Apply the fourth lifting step. */
3026                 lptr = &a[0];
3027                 hptr = &a[llen * stride];
3028                 if (parity) {
3029                         lptr2 = lptr;
3030                         hptr2 = hptr;
3031                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
3032                           ALPHA), lptr2[0]));
3033                         ++hptr2;
3034                         ++lptr2;
3035                         hptr += stride;
3036                 }
3037                 n = numrows - llen - parity - (parity == (numrows & 1));
3038                 while (n-- > 0) {
3039                         lptr2 = lptr;
3040                         hptr2 = hptr;
3041                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(ALPHA),
3042                           jpc_fix_add(lptr2[0], lptr2[stride])));
3043                         ++lptr2;
3044                         ++hptr2;
3045                         hptr += stride;
3046                         lptr += stride;
3047                 }
3048                 if (parity == (numrows & 1)) {
3049                         lptr2 = lptr;
3050                         hptr2 = hptr;
3051                         jpc_fix_minuseq(hptr2[0], jpc_fix_mul(jpc_dbltofix(2.0 *
3052                           ALPHA), lptr2[0]));
3053                         ++lptr2;
3054                         ++hptr2;
3055                 }
3056
3057         } else {
3058
3059 #if defined(WT_LENONE)
3060                 if (parity) {
3061                         lptr2 = &a[0];
3062                         lptr2[0] >>= 1;
3063                         ++lptr2;
3064                 }
3065 #endif
3066
3067         }
3068
3069 }
3070
3071 int jpc_ns_analyze(jpc_fix_t *a, int xstart, int ystart, int width, int height,
3072   int stride)
3073 {
3074
3075         int numrows = height;
3076         int numcols = width;
3077         int rowparity = ystart & 1;
3078         int colparity = xstart & 1;
3079         int i;
3080         jpc_fix_t *startptr;
3081         int maxcols;
3082
3083         maxcols = (numcols / JPC_QMFB_COLGRPSIZE) * JPC_QMFB_COLGRPSIZE;
3084         startptr = &a[0];
3085         for (i = 0; i < maxcols; i += JPC_QMFB_COLGRPSIZE) {
3086                 jpc_qmfb_split_colgrp(startptr, numrows, stride, rowparity);
3087                 jpc_ns_fwdlift_colgrp(startptr, numrows, stride, rowparity);
3088                 startptr += JPC_QMFB_COLGRPSIZE;
3089         }
3090         if (maxcols < numcols) {
3091                 jpc_qmfb_split_colres(startptr, numrows, numcols - maxcols, stride,
3092                   rowparity);
3093                 jpc_ns_fwdlift_colres(startptr, numrows, numcols - maxcols, stride,
3094                   rowparity);
3095         }
3096
3097         startptr = &a[0];
3098         for (i = 0; i < numrows; ++i) {
3099                 jpc_qmfb_split_row(startptr, numcols, colparity);
3100                 jpc_ns_fwdlift_row(startptr, numcols, colparity);
3101                 startptr += stride;
3102         }
3103
3104         return 0;
3105
3106 }
3107
3108 int jpc_ns_synthesize(jpc_fix_t *a, int xstart, int ystart, int width,
3109   int height, int stride)
3110 {
3111
3112         int numrows = height;
3113         int numcols = width;
3114         int rowparity = ystart & 1;
3115         int colparity = xstart & 1;
3116         int maxcols;
3117         jpc_fix_t *startptr;
3118         int i;
3119
3120         startptr = &a[0];
3121         for (i = 0; i < numrows; ++i) {
3122                 jpc_ns_invlift_row(startptr, numcols, colparity);
3123                 jpc_qmfb_join_row(startptr, numcols, colparity);
3124                 startptr += stride;
3125         }
3126
3127         maxcols = (numcols / JPC_QMFB_COLGRPSIZE) * JPC_QMFB_COLGRPSIZE;
3128         startptr = &a[0];
3129         for (i = 0; i < maxcols; i += JPC_QMFB_COLGRPSIZE) {
3130                 jpc_ns_invlift_colgrp(startptr, numrows, stride, rowparity);
3131                 jpc_qmfb_join_colgrp(startptr, numrows, stride, rowparity);
3132                 startptr += JPC_QMFB_COLGRPSIZE;
3133         }
3134         if (maxcols < numcols) {
3135                 jpc_ns_invlift_colres(startptr, numrows, numcols - maxcols, stride,
3136                   rowparity);
3137                 jpc_qmfb_join_colres(startptr, numrows, numcols - maxcols, stride,
3138                   rowparity);
3139         }
3140
3141         return 0;
3142
3143 }
3144