IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 7089


Ignore:
Timestamp:
May 8, 2006, 5:16:52 PM (20 years ago)
Author:
Paul Price
Message:

Added bounds checking to psImageSmooth. Discovered bugs: 1. The final
row and column were not being smoothed; and 2. The smoothed data was
being copied before the next smoothing, so that it could affect the
smoothing, resulting in a wider tail for the positive y axis. Fixed
bug 1 by removing the "-1"s. Fixed bug 2 by performing the copy after
the smoothing is performed, and carrying around an additional floating
row vector that cycles into the array of rows.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageConvolve.c

    r7077 r7089  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2006-05-05 23:55:56 $
     7 *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2006-05-09 03:16:52 $
    99 *
    1010 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    588588
    589589    // relevant terms
    590     int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side
     590    int Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
    591591    int Npixel = 2*Nrange + 1;          // Total number of pixels in convolution kernel
    592592    int Nx = image->numCols;            // Number of columns
    … …  
    596596case PS_TYPE_##TYPE: { \
    597597        /* generate normalized gaussian */ \
    598         psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_##TYPE); \
     598        psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \
    599599        gaussnorm->n = gaussnorm->nalloc; \
    600         double sum = 0.0; \
    601         double factor = -0.5/(sigma*sigma); \
    602         for (int i = -Nrange; i < Nrange + 1; i++) { \
    603             gaussnorm->data.TYPE[i+Nrange] = exp (factor*i*i); \
    604             sum += gaussnorm->data.TYPE[i+Nrange]; \
    605         } \
    606         for (int i = -Nrange; i < Nrange + 1; i++) { \
    607             gaussnorm->data.TYPE[i+Nrange] /= sum; \
     600        { \
     601            double sum = 0.0; \
     602            double factor = -0.5/(sigma*sigma); \
     603            for (int i = -Nrange; i < Nrange + 1; i++) { \
     604                gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \
     605                sum += gaussnorm->data.TYPE[i+Nrange]; \
     606            } \
     607            for (int i = -Nrange; i < Nrange + 1; i++) { \
     608                gaussnorm->data.TYPE[i+Nrange] /= sum; \
     609            } \
    608610        } \
    609611        ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \
    610612        \
    611         /* smooth in X direction */ \
    612         psVector *temp = psVectorAlloc (Nx, PS_TYPE_##TYPE); \
    613         temp->n = temp->nalloc; \
    614         for (int j = 0; j < Ny; j++) { \
    615             ps##TYPE *vi = image->data.TYPE[j]; \
    616             ps##TYPE *vo = temp->data.TYPE; \
    617             /* smooth first Nrange pixels, with renorm */ \
    618             int xMax = PS_MIN(Nrange, Nx); \
    619             for (int i = 0; i < xMax; i++, vi++, vo++) { \
    620                 ps##TYPE *vr = vi - i; \
    621                 ps##TYPE *vg = gauss - i; \
    622                 double g = 0; \
    623                 double s = 0; \
    624                 for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \
    625                     s += *vg * *vr; \
    626                     g += *vg; \
     613        /* Smooth in X direction */ \
     614        { \
     615            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
     616            calculation->n = calculation->nalloc; \
     617            for (int j = 0; j < Ny; j++) { \
     618                ps##TYPE *vi = image->data.TYPE[j]; \
     619                ps##TYPE *vo = calculation->data.TYPE; \
     620                int xMax = PS_MIN(Nrange, Nx); \
     621                int convRange = PS_MIN(Nrange + 1, Nx); \
     622                /* Smooth first Nrange pixels, with renorm */ \
     623                for (int i = 0; i < xMax; i++, vi++, vo++) { \
     624                    ps##TYPE *vr = vi - i; \
     625                    ps##TYPE *vg = gauss - i; \
     626                    double g = 0.0; \
     627                    double s = 0.0; \
     628                    for (int n = -i; n < convRange; n++, vr++, vg++) { \
     629                        s += *vg * *vr; \
     630                        g += *vg; \
     631                    } \
     632                    *vo = s / g; \
    627633                } \
    628                 *vo = s / g; \
    629             } \
    630             /* smooth middle pixels */ \
    631             for (int i = Nrange; i < Nx - Nrange; i++, vi++, vo++) { \
    632                 ps##TYPE *vr = vi - Nrange; \
    633                 ps##TYPE *vg = gauss - Nrange; \
    634                 double s = 0; \
    635                 for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \
    636                     s += *vg * *vr; \
     634                /* If that's all the pixels we have, then we're done already */ \
     635                if (Nx > Nrange) { \
     636                    /* Smooth middle pixels */ \
     637                    for (int i = Nrange; i < Nx - Nrange; i++, vi++, vo++) { \
     638                        ps##TYPE *vr = vi - Nrange; \
     639                        ps##TYPE *vg = gauss - Nrange; \
     640                        double s = 0; \
     641                        for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \
     642                            s += *vg * *vr; \
     643                        } \
     644                        *vo = s; \
     645                    } \
     646                    /* Smooth last Nrange pixels, with renorm */ \
     647                    /* XXX does this miss the last column? */ \
     648                    for (int i = Nx - Nrange; i < Nx; i++, vi++, vo++) { \
     649                        ps##TYPE *vr = vi - Nrange; \
     650                        ps##TYPE *vg = gauss - Nrange; \
     651                        double g = 0.0; \
     652                        double s = 0.0; \
     653                        for (int n = -Nrange; n < Nx - i; n++, vr++, vg++) { \
     654                            s += *vg * *vr; \
     655                            g += *vg; \
     656                        } \
     657                        *vo = s / g; \
     658                    } \
    637659                } \
    638                 *vo = s; \
    639             } \
    640             /* smooth last Nrange pixels, with renorm */ \
    641             /* XXX does this miss the last column? */ \
    642             for (int i = PS_MAX(Nx - Nrange, 0); i < Nx; i++, vi++, vo++) { \
    643                 ps##TYPE *vr = vi - Nrange; \
    644                 ps##TYPE *vg = gauss - Nrange; \
    645                 double g = 0; \
    646                 double s = 0; \
    647                 for (int n = -Nrange; n < Nx - i - 1; n++, vr++, vg++) { \
    648                     s += *vg * *vr; \
    649                     g += *vg; \
    650                 } \
    651                 *vo = s / g; \
    652             } \
    653             memcpy (image->data.TYPE[j], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
    654         } \
    655         psFree (temp); \
    656         \
    657         /* smooth in Y direction */ \
    658         psArray *temprows = psArrayAlloc (Nrange); \
    659         temprows->n = Nrange; \
     660                memcpy(image->data.TYPE[j], calculation->data.TYPE, Nx*sizeof(ps##TYPE)); \
     661            } \
     662            psFree(calculation); \
     663        } \
     664        \
     665        /* Smooth in Y direction */ \
     666        psArray *rows = psArrayAlloc(Nrange); \
     667        rows->n = Nrange; \
     668        /* Smooth the first Nrange pixels, with renorm */ \
    660669        int yMax = PS_MIN(Nrange, Ny); \
     670        int convRange = PS_MIN(Nrange + 1, Ny); \
    661671        for (int j = 0; j < yMax; j++) { \
    662             temp = psVectorAlloc (Nx, PS_TYPE_##TYPE); \
    663             temp->n = temp->nalloc; \
    664             /* zero the output row */ \
    665             memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
    666             for (int n = -j, sum = 0; n < Nrange + 1; n++) { sum += gauss[n]; } \
    667             for (int n = -j; n < Nrange + 1; n++) { \
     672            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
     673            calculation->n = calculation->nalloc; \
     674            /* Zero the output row */ \
     675            memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
     676            double sum = 0.0; \
     677            for (int n = -j; n < convRange; n++) { sum += gauss[n]; } \
     678            for (int n = -j; n < convRange; n++) { \
    668679                ps##TYPE *vi = image->data.TYPE[j+n]; \
    669                 ps##TYPE *vo = temp->data.TYPE; \
    670                 double g = gauss[n]; \
     680                ps##TYPE *vo = calculation->data.TYPE; \
     681                double g = gauss[n] / sum; \
    671682                for (int i = 0; i < Nx; i++, vi++, vo++) { \
    672683                    *vo += *vi * g; \
    673684                } \
    674685            } \
    675             /* save output rows on temp array */ \
    676             temprows->data[j] = temp; \
    677         } \
    678         for (int j = Nrange; j < Ny - Nrange; j++) { \
    679             /* save the Nrange-offset output row, then zero */ \
    680             int Nr = j % Nrange; \
    681             temp = temprows->data[Nr]; \
    682             memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
    683             memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
    684             for (int n = -Nrange; n < Nrange + 1; n++) { \
    685                 ps##TYPE *vi = image->data.TYPE[j+n]; \
    686                 ps##TYPE *vo = temp->data.TYPE; \
    687                 double g = gauss[n]; \
    688                 for (int i = 0; i < Nx; i++, vi++, vo++) { \
    689                     *vo += *vi * g; \
     686            /* Save output rows on temp array of rows */ \
     687            rows->data[j] = calculation; \
     688        } \
     689        if (Ny < Nrange) { \
     690            /* Need to save the first bit, then we're done */ \
     691            for (int j = 0; j < Ny; j++) {  \
     692                psVector *save = rows->data[j]; \
     693                memcpy(image->data.TYPE[j], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
     694            } \
     695        } else { \
     696            /* Smooth middle pixels */ \
     697            psVector *calculation = psVectorAlloc(Nx, PS_TYPE_##TYPE); \
     698            for (int j = Nrange; j < Ny - Nrange; j++) { \
     699                memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
     700                for (int n = -Nrange; n < Nrange + 1; n++) { \
     701                    ps##TYPE *vi = image->data.TYPE[j+n]; \
     702                    ps##TYPE *vo = calculation->data.TYPE; \
     703                    double g = gauss[n]; \
     704                    for (int i = 0; i < Nx; i++, vi++, vo++) { \
     705                        *vo += *vi * g; \
     706                    } \
    690707                } \
    691             } \
    692         } \
    693         for (int j = PS_MAX(Ny - Nrange, 0); j < Ny; j++) { \
    694             /* save the Nrange-offset output row, then zero */ \
    695             int Nr = j % Nrange; \
    696             temp = temprows->data[Nr]; \
    697             memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
    698             memset (temp->data.TYPE, 0, Nx*sizeof(PS_TYPE_##TYPE)); \
    699             for (int n = -Nrange, sum = 0; n < Ny - j - 1; n++) { sum += gauss[n]; } \
    700             for (int n = -Nrange; n < Ny - j - 1; n++) { \
    701                 ps##TYPE *vi = image->data.TYPE[j+n]; \
    702                 ps##TYPE *vo = temp->data.TYPE; \
    703                 double g = gauss[n]; \
    704                 for (int i = 0; i < Nx; i++, vi++, vo++) { \
    705                     *vo += *vi * g; \
     708                /* Write the output row */ \
     709                int Nr = j % Nrange; \
     710                psVector *save = rows->data[Nr]; \
     711                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
     712                /* Juggle the pointers, so that next run we use the one we just wrote */ \
     713                rows->data[Nr] = calculation; \
     714                calculation = save; \
     715            } \
     716            /* Smooth last Nrange pixels, with renorm */ \
     717            for (int j = Ny - Nrange; j < Ny; j++) { \
     718                /* save the Nrange-offset output row, then zero */ \
     719                memset(calculation->data.TYPE, 0, Nx*sizeof(ps##TYPE)); \
     720                double sum = 0.0; \
     721                for (int n = -Nrange; n < Ny - j; n++) { sum += gauss[n]; } \
     722                for (int n = -Nrange; n < Ny - j; n++) { \
     723                    ps##TYPE *vi = image->data.TYPE[j+n]; \
     724                    ps##TYPE *vo = calculation->data.TYPE; \
     725                    double g = gauss[n] / sum; \
     726                    for (int i = 0; i < Nx; i++, vi++, vo++) { \
     727                        *vo += *vi * g; \
     728                    } \
    706729                } \
    707             } \
    708         } \
    709         for (int j = Ny; j < Ny + Nrange; j++) { \
    710             /* save the Nrange-offset output row, then zero */ \
    711             int Nr = j % Nrange; \
    712             temp = temprows->data[Nr]; \
    713             memcpy (image->data.TYPE[j-Nrange], temp->data.TYPE, Nx*sizeof(ps##TYPE)); \
    714         } \
    715         psFree (temprows); \
    716         psFree (gaussnorm); \
     730                /* Write the output row */ \
     731                int Nr = j % Nrange; \
     732                psVector *save = rows->data[Nr]; \
     733                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
     734                /* Juggle the pointers, so that next run we use the one we just wrote */ \
     735                rows->data[Nr] = calculation; \
     736                calculation = save; \
     737            } \
     738            psFree(calculation); \
     739            /* Write the remaining rows */ \
     740            for (int j = Ny; j < Ny + Nrange; j++) {  \
     741                int Nr = j % Nrange; \
     742                psVector *save = rows->data[Nr]; \
     743                memcpy(image->data.TYPE[j-Nrange], save->data.TYPE, Nx*sizeof(ps##TYPE)); \
     744            } \
     745        } \
     746        psFree(rows); \
     747        psFree(gaussnorm); \
    717748        break; \
    718749    }
Note: See TracChangeset for help on using the changeset viewer.