wavelet.c

Go to the documentation of this file.
00001 /*                       W A V E L E T . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 2004-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @addtogroup wavelet */
00021 /** @{ */
00022 /** @file libbn/wavelet.c
00023  *
00024  * @brief
00025  *  This is a standard wavelet library that takes a given data buffer of some data
00026  *  type and then performs a wavelet transform on that data.
00027  *
00028  * The transform
00029  *  operations available are to either decompose or reconstruct a signal into its
00030  *  corresponding wavelet form based on the haar wavelet.
00031  *
00032  *  Wavelet decompose/reconstruct operations
00033  *
00034  *      - bn_wlt_haar_1d_double_decompose(tbuffer, buffer, dimen, channels, limit)
00035  *      - bn_wlt_haar_1d_float_decompose (tbuffer, buffer, dimen, channels, limit)
00036  *      - bn_wlt_haar_1d_char_decompose  (tbuffer, buffer, dimen, channels, limit)
00037  *      - bn_wlt_haar_1d_short_decompose (tbuffer, buffer, dimen, channels, limit)
00038  *      - bn_wlt_haar_1d_int_decompose   (tbuffer, buffer, dimen, channels, limit)
00039  *      - bn_wlt_haar_1d_long_decompose  (tbuffer, buffer, dimen, channels, limit)
00040  *
00041  *      - bn_wlt_haar_1d_double_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
00042  *      - bn_wlt_haar_1d_float_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
00043  *      - bn_wlt_haar_1d_char_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
00044  *      - bn_wlt_haar_1d_short_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
00045  *      - bn_wlt_haar_1d_int_reconstruct   (tbuffer, buffer, dimen, channels, sub_sz, limit)
00046  *      - bn_wlt_haar_1d_long_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
00047  *
00048  *      - bn_wlt_haar_2d_double_decompose(tbuffer, buffer, dimen, channels, limit)
00049  *      - bn_wlt_haar_2d_float_decompose (tbuffer, buffer, dimen, channels, limit)
00050  *      - bn_wlt_haar_2d_char_decompose  (tbuffer, buffer, dimen, channels, limit)
00051  *      - bn_wlt_haar_2d_short_decompose (tbuffer, buffer, dimen, channels, limit)
00052  *      - bn_wlt_haar_2d_int_decompose   (tbuffer, buffer, dimen, channels, limit)
00053  *      - bn_wlt_haar_2d_long_decompose  (tbuffer, buffer, dimen, channels, limit)
00054  *
00055  *      - bn_wlt_haar_2d_double_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
00056  *      - bn_wlt_haar_2d_float_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
00057  *      - bn_wlt_haar_2d_char_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
00058  *      - bn_wlt_haar_2d_short_reconstruct (tbuffer, buffer, dimen, channels, sub_sz, limit)
00059  *      - bn_wlt_haar_2d_int_reconstruct   (tbuffer, buffer, dimen, channels, sub_sz, limit)
00060  *      - bn_wlt_haar_2d_long_reconstruct  (tbuffer, buffer, dimen, channels, sub_sz, limit)
00061  *
00062  *      - bn_wlt_haar_2d_double_decompose2(tbuffer, buffer, width, height, channels, limit)
00063  *      - bn_wlt_haar_2d_float_decompose2 (tbuffer, buffer, width, height, channels, limit)
00064  *      - bn_wlt_haar_2d_char_decompose2  (tbuffer, buffer, width, height, channels, limit)
00065  *      - bn_wlt_haar_2d_short_decompose2 (tbuffer, buffer, width, height, channels, limit)
00066  *      - bn_wlt_haar_2d_int_decompose2   (tbuffer, buffer, width, height, channels, limit)
00067  *      - bn_wlt_haar_2d_long_decompose2  (tbuffer, buffer, width, height, channels, limit)
00068  *
00069  *      - bn_wlt_haar_2d_double_reconstruct2(tbuffer, buffer, width, height, channels, sub_sz, limit)
00070  *      - bn_wlt_haar_2d_float_reconstruct2 (tbuffer, buffer, width, height, channels, sub_sz, limit)
00071  *      - bn_wlt_haar_2d_char_reconstruct2  (tbuffer, buffer, width, height, channels, sub_sz, limit)
00072  *      - bn_wlt_haar_2d_short_reconstruct2 (tbuffer, buffer, width, height, channels, sub_sz, limit)
00073  *      - bn_wlt_haar_2d_int_reconstruct2   (tbuffer, buffer, width, height, channels, sub_sz, limit)
00074  *      - bn_wlt_haar_2d_long_reconstruct2  (tbuffer, buffer, width, height, channels, sub_sz, limit)
00075  *
00076  *
00077  *  For greatest accuracy, it is preferable to convert everything to "double"
00078  *  and decompose/reconstruct with that.  However, there are useful
00079  *  properties to performing the decomposition and/or reconstruction in
00080  *  various data types (most notably char).
00081  *
00082  *  Rather than define all of these routines explicitly, we define
00083  *  2 macros "decompose" and "reconstruct" which embody the structure of
00084  *  the function (which is common to all of them).  We then instatiate
00085  *  these macros once for each of the data types.  It's ugly, but it
00086  *  assures that a change to the structure of one operation type
00087  *  (decompose or reconstruct) occurs for all data types.
00088  *
00089  *
00090  *
00091  *
00092  *  bn_wlt_haar_1d_*_decompose(tbuffer, buffer, dimen, channels, limit)
00093  *  @par Parameters:
00094  *      - tbuffer     a temporary data buffer 1/2 as large as "buffer". See (1) below.
00095  *      - buffer      pointer to the data to be decomposed
00096  *      - dimen    the number of samples in the data buffer
00097  *      - channels the number of values per sample
00098  *      - limit    the extent of the decomposition
00099  *
00100  *  Perform a Haar wavelet decomposition on the data in buffer "buffer".  The
00101  *  decomposition is done "in place" on the data, hence the values in "buffer"
00102  *  are not preserved, but rather replaced by their decomposition.
00103  *  The number of original samples in the buffer (parameter "dimen") and the
00104  *  decomposition limit ("limit") must both be a power of 2 (e.g. 512, 1024).
00105  *  The buffer is decomposed into "average" and "detail" halves until the
00106  *  size of the "average" portion reaches "limit".  Simultaneous
00107  *  decomposition of multi-plane (e.g. pixel) data, can be performed by
00108  *  indicating the number of planes in the "channels" parameter.
00109  *
00110  *  (1) The process requires a temporary buffer which is 1/2 the size of the
00111  *  longest span to be decomposed.  If the "tbuffer" argument is non-null then
00112  *  it is a pointer to a temporary buffer.  If the pointer is NULL, then a
00113  *  local temporary buffer will be allocated (and freed).
00114  *
00115  *  @par Examples:
00116  @code
00117  double dbuffer[512], cbuffer[256];
00118  ...
00119  bn_wlt_haar_1d_double_decompose(cbuffer, dbuffer, 512, 1, 1);
00120  @endcode
00121  *
00122  *    performs complete decomposition on the data in array "dbuffer".
00123  *
00124  @code
00125  double buffer[3][512];  /_* 512 samples, 3 values/sample (e.g. RGB?)*_/
00126  double tbuffer[3][256];         /_* the temporary buffer *_/
00127  ...
00128  bn_wlt_haar_1d_double_decompose(tbuffer, buffer, 512, 3, 1);
00129  @endcode
00130  *
00131  *    This will completely decompose the data in buffer.  The first sample will
00132  *    be the average of all the samples.  Alternatively:
00133  *
00134  * bn_wlt_haar_1d_double_decompose(tbuffer, buffer, 512, 3, 64);
00135  *
00136  *    decomposes buffer into a 64-sample "average image" and 3 "detail" sets.
00137  *
00138  *  bn_wlt_haar_1d_*_reconstruct(tbuffer, buffer, dimen, channels, sub_sz, limit)
00139  *
00140  */
00141 /** @} */
00142 
00143 #include "common.h"
00144 
00145 #include <stdio.h>
00146 #include <string.h>
00147 
00148 #include "bu.h"
00149 #include "vmath.h"
00150 #include "bn.h"
00151 
00152 /**
00153  * This source file uses C-styple "templates" where functions specific
00154  * to a set of specified data types are automatically
00155  * declared/provided via one single macro implementation.
00156  */
00157 #define decompose_1d(DATATYPE) bn_wlt_haar_1d_ ## DATATYPE ## _decompose
00158 
00159 
00160 /**
00161  * DATATYPE *tbuffer;      // temporary buffer
00162  * DATATYPE *buffer;       // data buffer
00163  * unsigned long dimen;    // # of samples in data buffer
00164  * unsigned long channels; // # of data values per sample
00165  * unsigned long limit;    // extent of decomposition
00166  */
00167 
00168 #define make_wlt_haar_1d_decompose(DATATYPE)                            \
00169     void                                                                \
00170     decompose_1d(DATATYPE)                                              \
00171         ( DATATYPE *tbuffer, DATATYPE *buffer, unsigned long dimen, unsigned long channels, unsigned long limit ) \
00172     {                                                                   \
00173         register DATATYPE *detail;                                      \
00174         register DATATYPE *avg;                                         \
00175         unsigned long img_size;                                         \
00176         unsigned long half_size;                                        \
00177         int do_free = 0;                                                \
00178         unsigned long x, x_tmp, d, i, j_idx;                            \
00179         register fastf_t onehalf = (fastf_t)0.5;                        \
00180                                                                         \
00181         CK_POW_2( dimen );                                              \
00182                                                                         \
00183         if ( ! tbuffer ) {                                              \
00184             tbuffer = (DATATYPE *)bu_malloc(                            \
00185                 (dimen/2) * channels * sizeof( *buffer ),               \
00186                 "1d wavelet buffer");                                   \
00187             do_free = 1;                                                \
00188         }                                                               \
00189                                                                         \
00190         /* each iteration of this loop decomposes the data into 2 halves: \
00191          * the "average image" and the "image detail"                   \
00192          */                                                             \
00193         for (img_size = dimen; img_size > limit; img_size = half_size ) { \
00194                                                                         \
00195             half_size = img_size/2;                                     \
00196                                                                         \
00197             detail = tbuffer;                                           \
00198             avg = buffer;                                               \
00199                                                                         \
00200             for ( x=0; x < img_size; x += 2 ) {                         \
00201                 x_tmp = x*channels;                                     \
00202                                                                         \
00203                 for (d=0; d < channels; d++, avg++, detail++) {         \
00204                     i = x_tmp + d;                                      \
00205                     j_idx = i + channels;                               \
00206                     *detail = (buffer[i] - buffer[j_idx]) * onehalf;    \
00207                     *avg    = (buffer[i] + buffer[j_idx]) * onehalf;    \
00208                 }                                                       \
00209             }                                                           \
00210                                                                         \
00211             /* "avg" now points to the first element AFTER the "average \
00212              * image" section, and hence is the START of the "image     \
00213              * detail" portion.  Convenient, since we now want to copy  \
00214              * the contents of "tbuffer" (which holds the image detail) into \
00215              * place.                                                   \
00216              */                                                         \
00217             memcpy(avg, tbuffer, sizeof(*buffer) * channels * half_size); \
00218         }                                                               \
00219                                                                         \
00220         if (do_free)                                                    \
00221             bu_free( (genptr_t)tbuffer, "1d wavelet buffer");           \
00222     }
00223 
00224 
00225 #define reconstruct(DATATYPE ) bn_wlt_haar_1d_ ## DATATYPE ## _reconstruct
00226 
00227 #define make_wlt_haar_1d_reconstruct( DATATYPE )                        \
00228     void                                                                \
00229     reconstruct(DATATYPE)                                               \
00230         ( DATATYPE *tbuffer, DATATYPE *buffer, unsigned long dimen, unsigned long channels, unsigned long subimage_size, unsigned long limit ) \
00231     {                                                                   \
00232         register DATATYPE *detail;                                      \
00233         register DATATYPE *avg;                                         \
00234         unsigned long img_size;                                         \
00235         unsigned long dbl_size;                                         \
00236         int do_free = 0;                                                \
00237         unsigned long x_tmp, d, x, i, j_idx;                            \
00238                                                                         \
00239         CK_POW_2( subimage_size );                                      \
00240         CK_POW_2( dimen );                                              \
00241         CK_POW_2( limit );                                              \
00242                                                                         \
00243         if ( ! (subimage_size < dimen) ) {                              \
00244             bu_log("%s:%d Dimension %lu should be greater than subimage size (%lu)\n", \
00245                    __FILE__, __LINE__, dimen, subimage_size);           \
00246             bu_bomb("reconstruct");                                     \
00247         }                                                               \
00248                                                                         \
00249         if ( ! (subimage_size < limit) ) {                              \
00250             bu_log("%s:%d Channels limit %lu should be greater than subimage size (%lu)\n", \
00251                    __FILE__, __LINE__, limit, subimage_size);           \
00252             bu_bomb("reconstruct");                                     \
00253         }                                                               \
00254                                                                         \
00255         if ( ! (limit <= dimen) ) {                                     \
00256             bu_log("%s:%d Dimension %lu should be greater than or equal to the channels limit (%lu)\n", \
00257                    __FILE__, __LINE__, dimen, limit);                   \
00258             bu_bomb("reconstruct");                                     \
00259         }                                                               \
00260                                                                         \
00261                                                                         \
00262         if ( ! tbuffer ) {                                              \
00263             tbuffer = ( DATATYPE *)bu_malloc((dimen/2) * channels * sizeof( *buffer ), \
00264                                              "1d wavelet reconstruct tmp buffer"); \
00265             do_free = 1;                                                \
00266         }                                                               \
00267                                                                         \
00268         /* Each iteration of this loop reconstructs an image twice as   \
00269          * large as the original using a "detail image".                \
00270          */                                                             \
00271         for (img_size=subimage_size; img_size < limit; img_size=dbl_size) { \
00272             dbl_size = img_size * 2;                                    \
00273                                                                         \
00274             d = img_size * channels;                                    \
00275             detail = &buffer[ d ];                                      \
00276                                                                         \
00277             /* copy the original or "average" data to temporary buffer */ \
00278             avg = tbuffer;                                              \
00279             memcpy(avg, buffer, sizeof(*buffer) * d);                   \
00280                                                                         \
00281                                                                         \
00282             for (x=0; x < dbl_size; x += 2 ) {                          \
00283                 x_tmp = x * channels;                                   \
00284                 for (d=0; d < channels; d++, avg++, detail++ ) {        \
00285                     i = x_tmp + d;                                      \
00286                     j_idx = i + channels;                               \
00287                     buffer[i] = *avg + *detail;                         \
00288                     buffer[j_idx] = *avg - *detail;                     \
00289                 }                                                       \
00290             }                                                           \
00291         }                                                               \
00292                                                                         \
00293         if (do_free)                                                    \
00294             bu_free( (genptr_t)tbuffer,                                 \
00295                      "1d wavelet reconstruct tmp buffer");              \
00296     }
00297 
00298 /* Believe it or not, this is where the actual code is generated */
00299 
00300 make_wlt_haar_1d_decompose(double)
00301 make_wlt_haar_1d_reconstruct(double)
00302 
00303 make_wlt_haar_1d_decompose(float)
00304 make_wlt_haar_1d_reconstruct(float)
00305 
00306 make_wlt_haar_1d_decompose(char)
00307 make_wlt_haar_1d_reconstruct(char)
00308 
00309 make_wlt_haar_1d_decompose(int)
00310 make_wlt_haar_1d_reconstruct(int)
00311 
00312 make_wlt_haar_1d_decompose(short)
00313 make_wlt_haar_1d_reconstruct(short)
00314 
00315 make_wlt_haar_1d_decompose(long)
00316 make_wlt_haar_1d_reconstruct(long)
00317 
00318 
00319 #define decompose_2d( DATATYPE ) bn_wlt_haar_2d_ ## DATATYPE ## _decompose
00320 
00321 #define make_wlt_haar_2d_decompose(DATATYPE)                            \
00322     void                                                                \
00323     decompose_2d(DATATYPE)                                              \
00324         (DATATYPE *tbuffer, DATATYPE *buffer, unsigned long dimen, unsigned long channels, unsigned long limit) \
00325     {                                                                   \
00326         register DATATYPE *detail;                                      \
00327         register DATATYPE *avg;                                         \
00328         unsigned long img_size;                                         \
00329         unsigned long half_size;                                        \
00330         unsigned long x, y, x_tmp, y_tmp, d, i, j_idx;                  \
00331         register fastf_t onehalf = (fastf_t)0.5;                        \
00332                                                                         \
00333         CK_POW_2( dimen );                                              \
00334                                                                         \
00335         if ( ! tbuffer ) {                                              \
00336             tbuffer = (DATATYPE *)bu_malloc(                            \
00337                 (dimen/2) * channels * sizeof( *buffer ),               \
00338                 "1d wavelet buffer");                                   \
00339         }                                                               \
00340                                                                         \
00341         /* each iteration of this loop decomposes the data into 4 quarters: \
00342          * the "average image", the horizontal detail, the vertical detail \
00343          * and the horizontal-vertical detail                           \
00344          */                                                             \
00345         for (img_size = dimen; img_size > limit; img_size = half_size ) { \
00346             half_size = img_size/2;                                     \
00347                                                                         \
00348             /* do a horizontal detail decomposition first */            \
00349             for (y=0; y < img_size; y++ ) {                             \
00350                 y_tmp = y * dimen * channels;                           \
00351                                                                         \
00352                 detail = tbuffer;                                       \
00353                 avg = &buffer[y_tmp];                                   \
00354                                                                         \
00355                 for (x=0; x < img_size; x += 2 ) {                      \
00356                     x_tmp = x*channels + y_tmp;                         \
00357                                                                         \
00358                     for (d=0; d < channels; d++, avg++, detail++) {     \
00359                         i = x_tmp + d;                                  \
00360                         j_idx = i + channels;                           \
00361                         *detail = (buffer[i] - buffer[j_idx]) * onehalf; \
00362                         *avg    = (buffer[i] + buffer[j_idx]) * onehalf; \
00363                     }                                                   \
00364                 }                                                       \
00365                 /* "avg" now points to the first element AFTER the      \
00366                  * "average image" section, and hence is the START      \
00367                  * of the "image detail" portion.  Convenient, since    \
00368                  * we now want to copy the contents of "tbuffer" (which \
00369                  * holds the image detail) into place.                  \
00370                  */                                                     \
00371                 memcpy(avg, tbuffer, sizeof(*buffer) * channels * half_size); \
00372             }                                                           \
00373                                                                         \
00374             /* Now do the vertical decomposition */                     \
00375             for (x=0; x < img_size; x ++ ) {                            \
00376                 x_tmp = x*channels;                                     \
00377                                                                         \
00378                 detail = tbuffer;                                       \
00379                 avg = &buffer[x_tmp];                                   \
00380                                                                         \
00381                 for (y=0; y < img_size; y += 2) {                       \
00382                     y_tmp =y*dimen*channels + x_tmp;                    \
00383                                                                         \
00384                     for (d=0; d < channels; d++, avg++, detail++) {     \
00385                         i = y_tmp + d;                                  \
00386                         j_idx = i + dimen*channels;                     \
00387                         *detail = (buffer[i] - buffer[j_idx]) * onehalf; \
00388                         *avg    = (buffer[i] + buffer[j_idx]) * onehalf; \
00389                     }                                                   \
00390                     avg += (dimen-1)*channels;                          \
00391                 }                                                       \
00392                                                                         \
00393                 /* "avg" now points to the element ABOVE the            \
00394                  * last "average image" pixel or the first "detail"     \
00395                  * location in the user buffer.                         \
00396                  *                                                      \
00397                  * There is no memcpy for the columns, so we have to    \
00398                  * copy the data back to the user buffer ourselves.     \
00399                  */                                                     \
00400                 detail = tbuffer;                                       \
00401                 for (y=half_size; y < img_size; y++) {                  \
00402                     for (d=0; d < channels; d++) {                      \
00403                         *avg++ = *detail++;                             \
00404                     }                                                   \
00405                     avg += (dimen-1)*channels;                          \
00406                 }                                                       \
00407             }                                                           \
00408         }                                                               \
00409     }
00410 
00411 
00412 #define reconstruct_2d( DATATYPE ) bn_wlt_haar_2d_ ## DATATYPE ## _reconstruct
00413 
00414 #define make_wlt_haar_2d_reconstruct(DATATYPE)                          \
00415     void                                                                \
00416     reconstruct_2d(DATATYPE)                                            \
00417         (DATATYPE *tbuf, DATATYPE *buf, unsigned long width, unsigned long channels, unsigned long avg_size, unsigned long limit) \
00418     {                                                                   \
00419         register DATATYPE *detail;                                      \
00420         register DATATYPE *avg;                                         \
00421         unsigned long img_size;                                         \
00422         unsigned long dbl_size;                                         \
00423         unsigned long x_tmp, d, x, i, j_idx;                            \
00424         unsigned long y, row_len, row_start;                            \
00425                                                                         \
00426         CK_POW_2( avg_size );                                           \
00427         CK_POW_2( width );                                              \
00428         CK_POW_2( limit );                                              \
00429                                                                         \
00430         /* XXX check for:                                               \
00431          * subimage_size < dimen && subimage_size < limit               \
00432          * limit <= dimen                                               \
00433          */                                                             \
00434                                                                         \
00435                                                                         \
00436         if ( ! tbuf ) {                                                 \
00437             tbuf = ( DATATYPE *)bu_malloc((width/2) * channels * sizeof( *buf ), \
00438                                           "1d wavelet reconstruct tmp buffer"); \
00439         }                                                               \
00440                                                                         \
00441         row_len = width * channels;                                     \
00442                                                                         \
00443         /* Each iteration of this loop reconstructs an image twice as   \
00444          * large as the original using a "detail image".                \
00445          */                                                             \
00446                                                                         \
00447         for (img_size = avg_size; img_size < limit; img_size = dbl_size) { \
00448             dbl_size = img_size * 2;                                    \
00449                                                                         \
00450                                                                         \
00451             /* first is a vertical reconstruction */                    \
00452             for (x=0; x < dbl_size; x++ ) {                             \
00453                 /* reconstruct column x */                              \
00454                                                                         \
00455                 /* copy column of "average" data to tbuf */             \
00456                 x_tmp = x*channels;                                     \
00457                 for (y=0; y < img_size; y++) {                          \
00458                     i = x_tmp + y*row_len;                              \
00459                     j_idx = y * channels;                               \
00460                     for (d=0; d < channels; d++) {                      \
00461                         tbuf[j_idx++] = buf[i++];                       \
00462                     }                                                   \
00463                 }                                                       \
00464                 avg = tbuf;                                             \
00465                 detail = &buf[x_tmp + img_size*row_len];                \
00466                                                                         \
00467                 /* reconstruct column */                                \
00468                 for (y=0; y < dbl_size; y += 2) {                       \
00469                                                                         \
00470                     i = x_tmp + y*row_len;                              \
00471                     j_idx = i + row_len;                                \
00472                                                                         \
00473                     for (d=0; d < channels; d++, avg++, detail++) {     \
00474                         buf[i++] = *avg + *detail;                      \
00475                         buf[j_idx++] = *avg - *detail;                  \
00476                     }                                                   \
00477                     detail += row_len - channels;                       \
00478                 }                                                       \
00479             }                                                           \
00480                                                                         \
00481             /* now a horizontal reconstruction */                       \
00482             for (y=0; y < dbl_size; y++ ) {                             \
00483                 /* reconstruct row y */                                 \
00484                                                                         \
00485                 /* copy "average" row to tbuf and set pointer to        \
00486                  * begining of "detail"                                 \
00487                  */                                                     \
00488                 d = img_size * channels;                                \
00489                 row_start = y*row_len;                                  \
00490                                                                         \
00491                                                                         \
00492                 avg = &buf[ row_start ];                                \
00493                 detail = &buf[ row_start + d];                          \
00494                                                                         \
00495                 memcpy(tbuf, avg, sizeof(*buf) * d);                    \
00496                 avg = tbuf;                                             \
00497                                                                         \
00498                 /* reconstruct row */                                   \
00499                 for (x=0; x < dbl_size; x += 2 ) {                      \
00500                     x_tmp = x * channels;                               \
00501                     i = row_start + x * channels;                       \
00502                     j_idx = i + channels;                               \
00503                                                                         \
00504                     for (d=0; d < channels; d++, avg++, detail++) {     \
00505                         buf[i++] = *avg + *detail;                      \
00506                         buf[j_idx++] = *avg - *detail;                  \
00507                     }                                                   \
00508                 }                                                       \
00509             }                                                           \
00510         }                                                               \
00511     }
00512 
00513 make_wlt_haar_2d_decompose(double)
00514 make_wlt_haar_2d_decompose(float)
00515 make_wlt_haar_2d_decompose(char)
00516 make_wlt_haar_2d_decompose(int)
00517 make_wlt_haar_2d_decompose(short)
00518 make_wlt_haar_2d_decompose(long)
00519 
00520 make_wlt_haar_2d_reconstruct(double)
00521 make_wlt_haar_2d_reconstruct(float)
00522 make_wlt_haar_2d_reconstruct(char)
00523 make_wlt_haar_2d_reconstruct(int)
00524 make_wlt_haar_2d_reconstruct(short)
00525 make_wlt_haar_2d_reconstruct(long)
00526 
00527 
00528 #define decompose_2d_2( DATATYPE ) bn_wlt_haar_2d_ ## DATATYPE ## _decompose2
00529 
00530 #define make_wlt_haar_2d_decompose2(DATATYPE)                           \
00531     void                                                                \
00532     decompose_2d_2(DATATYPE)                                            \
00533         (DATATYPE *tbuffer, DATATYPE *buffer, unsigned long width, unsigned long height, unsigned long channels, unsigned long limit) \
00534     {                                                                   \
00535         register DATATYPE *detail;                                      \
00536         register DATATYPE *avg;                                         \
00537         unsigned long img_wsize;                                        \
00538         unsigned long img_hsize;                                        \
00539         unsigned long half_wsize;                                       \
00540         unsigned long half_hsize;                                       \
00541         unsigned long x, y, x_tmp, y_tmp, d, i, j_idx;                  \
00542         register fastf_t onehalf = (fastf_t)0.5;                        \
00543                                                                         \
00544         CK_POW_2( width );                                              \
00545         CK_POW_2( height );                                             \
00546                                                                         \
00547         /* create a temp buffer the half the size of the larger dimension \
00548          */                                                             \
00549         if ( ! tbuffer ) {                                              \
00550             tbuffer = (DATATYPE *)bu_malloc(                            \
00551                 (((width>height)?width:height)/2) * channels * sizeof( *buffer ), \
00552                 "1d wavelet buffer");                                   \
00553         }                                                               \
00554                                                                         \
00555         /* each iteration of this loop decomposes the data into 4 quarters: \
00556          * the "average image", the horizontal detail, the vertical detail \
00557          * and the horizontal-vertical detail                           \
00558          */                                                             \
00559         for (img_wsize = width, img_hsize = height; (img_wsize > limit) && (img_hsize > limit); img_wsize = half_wsize, img_hsize = half_hsize ) { \
00560             half_wsize = img_wsize/2;                                   \
00561             half_hsize = img_hsize/2;                                   \
00562                                                                         \
00563             /* do a horizontal detail decomposition first */            \
00564             for (y=0; y < img_hsize; y++ ) {                            \
00565                 y_tmp = y * width * channels;                           \
00566                                                                         \
00567                 detail = tbuffer;                                       \
00568                 avg = &buffer[y_tmp];                                   \
00569                                                                         \
00570                 for (x=0; x < img_wsize; x += 2 ) {                     \
00571                     x_tmp = x*channels + y_tmp;                         \
00572                                                                         \
00573                     for (d=0; d < channels; d++, avg++, detail++) {     \
00574                         i = x_tmp + d;                                  \
00575                         j_idx = i + channels;                           \
00576                         *detail = (buffer[i] - buffer[j_idx]) * onehalf; \
00577                         *avg    = (buffer[i] + buffer[j_idx]) * onehalf; \
00578                     }                                                   \
00579                 }                                                       \
00580                 /* "avg" now points to the first element AFTER the      \
00581                  * "average image" section, and hence is the START      \
00582                  * of the "image detail" portion.  Convenient, since    \
00583                  * we now want to copy the contents of "tbuffer" (which \
00584                  * holds the image detail) into place.                  \
00585                  */                                                     \
00586                 memcpy(avg, tbuffer, sizeof(*buffer) * channels * half_wsize); \
00587             }                                                           \
00588                                                                         \
00589             /* Now do the vertical decomposition */                     \
00590             for (x=0; x < img_wsize; x ++ ) {                           \
00591                 x_tmp = x*channels;                                     \
00592                                                                         \
00593                 detail = tbuffer;                                       \
00594                 avg = &buffer[x_tmp];                                   \
00595                                                                         \
00596                 for (y=0; y < img_hsize; y += 2) {                      \
00597                     y_tmp =y*width*channels + x_tmp;                    \
00598                                                                         \
00599                     for (d=0; d < channels; d++, avg++, detail++) {     \
00600                         i = y_tmp + d;                                  \
00601                         j_idx = i + width*channels;                     \
00602                         *detail = (buffer[i] - buffer[j_idx]) * onehalf; \
00603                         *avg    = (buffer[i] + buffer[j_idx]) * onehalf; \
00604                     }                                                   \
00605                     avg += (width-1)*channels;                          \
00606                 }                                                       \
00607                                                                         \
00608                 /* "avg" now points to the element ABOVE the            \
00609                  * last "average image" pixel or the first "detail"     \
00610                  * location in the user buffer.                         \
00611                  *                                                      \
00612                  * There is no memcpy for the columns, so we have to    \
00613                  * copy the data back to the user buffer ourselves.     \
00614                  */                                                     \
00615                 detail = tbuffer;                                       \
00616                 for (y=half_hsize; y < img_hsize; y++) {                \
00617                     for (d=0; d < channels; d++) {                      \
00618                         *avg++ = *detail++;                             \
00619                     }                                                   \
00620                     avg += (width-1)*channels;                          \
00621                 }                                                       \
00622             }                                                           \
00623         }                                                               \
00624     }
00625 
00626 make_wlt_haar_2d_decompose2(double)
00627 make_wlt_haar_2d_decompose2(float)
00628 make_wlt_haar_2d_decompose2(char)
00629 make_wlt_haar_2d_decompose2(int)
00630 make_wlt_haar_2d_decompose2(short)
00631 make_wlt_haar_2d_decompose2(long)
00632 
00633 
00634 /*
00635  * Local Variables:
00636  * mode: C
00637  * tab-width: 8
00638  * indent-tabs-mode: t
00639  * c-file-style: "stroustrup"
00640  * End:
00641  * ex: shiftwidth=4 tabstop=8
00642  */
Generated on Tue Dec 11 13:14:28 2012 for LIBBN by  doxygen 1.6.3