IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36319


Ignore:
Timestamp:
Nov 28, 2013, 4:10:42 AM (13 years ago)
Author:
eugene
Message:

extended 2d cached-convolution to 15 pix radius

Location:
branches/eam_branches/ipp-20130904/psLib
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c

    r36307 r36319  
    2727
    2828// fixed set of radii and number of entries in circularly symmetric profile
    29 # define NRAD_MAX 44
    30 static float radii2[NRAD_MAX] = {  0.0,  1.0,  2.0,  4.0,  5.0,
    31                                    8.0,  9.0, 10.0, 13.0, 16.0,
    32                                    17.0, 18.0, 20.0, 25.0, 26.0,
    33                                    29.0, 32.0, 34.0, 36.0, 37.0,
    34                                    40.0, 41.0, 45.0, 49.0, 50.0,
    35                                    52.0, 53.0, 58.0, 61.0, 64.0,
    36                                    65.0, 68.0, 72.0, 73.0, 74.0,
    37                                    80.0, 81.0, 82.0, 85.0, 89.0,
    38                                    90.0, 97.0, 98.0, 100.0
     29# define NRAD_MAX 87
     30# define  RAD_MAX 15
     31static float radii2[NRAD_MAX] = {   0.0,   1.0,   2.0,   4.0,   5.0,
     32                                    8.0,   9.0,  10.0,  13.0,  16.0,
     33                                    17.0,  18.0,  20.0,  25.0,  26.0,
     34
     35                                    29.0,  32.0,  34.0,  36.0,  37.0,
     36                                    40.0,  41.0,  45.0,  49.0,  50.0,
     37                                    52.0,  53.0,  58.0,  61.0,  64.0,
     38
     39                                    65.0,  68.0,  72.0,  73.0,  74.0,
     40                                    80.0,  81.0,  82.0,  85.0,  89.0,
     41                                    90.0,  97.0,  98.0, 100.0, 101.0,
     42
     43                                    104.0, 106.0, 109.0, 113.0, 116.0,
     44                                    117.0, 121.0, 122.0, 125.0, 128.0,
     45                                    130.0, 136.0, 137.0, 144.0, 145.0,
     46
     47                                    146.0, 148.0, 149.0, 153.0, 157.0,
     48                                    160.0, 162.0, 164.0, 169.0, 170.0,
     49                                    173.0, 178.0, 180.0, 181.0, 185.0,
     50
     51                                    193.0, 194.0, 196.0, 197.0, 200.0,
     52                                    202.0, 205.0, 208.0, 212.0, 218.0,
     53                                    221.0, 225.0
    3954};
    4055
    … …  
    4257                                     4,  4, 8,  8,  4,
    4358                                     8,  4, 8, 12,  8,
     59
    4460                                     8,  4, 8,  4,  8,
    4561                                     8,  8, 8,  4, 12,
    4662                                     8,  8, 8,  8,  4,
     63
    4764                                     16, 8, 4,  8,  8,
    48                                       8, 4, 8, 16,  8,
    49                                       8, 8, 4, 12
     65                                     8, 4, 8, 16,  8,
     66                                     8, 8, 4, 12, 8,
     67
     68                                     8, 8, 8, 8, 8,
     69                                     8, 4, 8, 16, 4,
     70                                     16, 8, 8, 8, 16,
     71
     72                                     8, 8, 8, 8, 8,
     73                                     8, 4, 8, 12, 16,
     74                                     8, 8, 8, 8, 16,
     75
     76                                     8, 8, 4, 8, 12,
     77                                     8, 16, 8, 8, 8,
     78                                     16, 12
    5079};
    5180
    … …  
    85114psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma) {
    86115
     116    psAssert (isfinite(Nsigma) && (Nsigma > 0.5), "Nsigma is not valid");
     117
    87118    psImageSmooth2dCacheData *smdata = psAlloc(sizeof(psImageSmooth2dCacheData));
    88119    psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth2dCacheDataFree);
    … …  
    100131    // check for NULL structure elements?
    101132    int Ns = (int)(smdata->Nsigma * sigma);
    102     Ns = PS_MAX (3, PS_MIN (Ns, 10));
     133    Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX));
    103134    smdata->Ns = Ns;
    104135
    … …  
    108139    smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX);
    109140
    110     float sum = 0.0;
     141    double sum = 0.0;
    111142    for (int i = 0; i < NRAD_MAX; i++) {
    112143        if (radii2[i] > Ns2) {
    … …  
    130161    // check for NULL structure elements?
    131162    int Ns = (int)(smdata->Nsigma * sigma);
    132     Ns = PS_MAX (3, PS_MIN (Ns, 10));
     163    Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX));
    133164    smdata->Ns = Ns;
    134165
    … …  
    275306            if (Ns <= 10) goto finish;
    276307
     308            ADD_RAND ( 44,  1, 10);  // r^2 = 101
     309            ADD_RAND ( 45,  2, 10);  // r^2 = 104
     310            ADD_RAND ( 46,  5,  9);  // r^2 = 106
     311            ADD_RAND ( 47,  3, 10);  // r^2 = 109
     312            ADD_RAND ( 48,  7,  8);  // r^2 = 113
     313            ADD_RAND ( 49,  4, 10);  // r^2 = 116
     314            ADD_RAND ( 50,  6,  9);  // r^2 = 117
     315            ADD_AXIS ( 51, 11)    ;  // r^2 = 121
     316            if (Ns <= 11) goto finish;
     317
     318            ADD_RAND ( 52,  1, 11);  // r^2 = 122
     319            ADD_RAND ( 53,  2, 11);  // r^2 = 125 *
     320            ADD_RAND ( 53,  5, 10);  // r^2 = 125 *
     321            ADD_DIAG ( 54,  8)    ;  // r^2 = 128
     322            ADD_RAND ( 55,  3, 11);  // r^2 = 130 *
     323            ADD_RAND ( 55,  7,  9);  // r^2 = 130 *
     324            ADD_RAND ( 56,  6, 10);  // r^2 = 136
     325            ADD_RAND ( 57,  4, 11);  // r^2 = 137
     326            ADD_AXIS ( 58, 12)    ;  // r^2 = 144
     327            if (Ns <= 12) goto finish;
     328
     329            ADD_RAND ( 59,  1, 12);  // r^2 = 145 *
     330            ADD_RAND ( 59,  8,  9);  // r^2 = 145 *
     331            ADD_RAND ( 60,  5, 11);  // r^2 = 146
     332            ADD_RAND ( 61,  2, 12);  // r^2 = 148
     333            ADD_RAND ( 62,  7, 10);  // r^2 = 149
     334            ADD_RAND ( 63,  3, 12);  // r^2 = 153
     335            ADD_RAND ( 64,  6, 11);  // r^2 = 157
     336            ADD_RAND ( 65,  4, 12);  // r^2 = 160
     337            ADD_DIAG ( 66,  9)    ;  // r^2 = 162
     338            ADD_RAND ( 67,  8, 10);  // r^2 = 164
     339            ADD_RAND ( 68,  5, 12);  // r^2 = 169 *
     340            ADD_AXIS ( 68, 13)    ;  // r^2 = 169 *
     341            if (Ns <= 13) goto finish;
     342
     343            ADD_RAND ( 69,  1, 13);  // r^2 = 170 *
     344            ADD_RAND ( 69,  7, 11);  // r^2 = 170 *
     345            ADD_RAND ( 70,  2, 13);  // r^2 = 173
     346            ADD_RAND ( 71,  3, 13);  // r^2 = 178
     347            ADD_RAND ( 72,  6, 12);  // r^2 = 180
     348            ADD_RAND ( 73,  9, 10);  // r^2 = 181
     349            ADD_RAND ( 74,  4, 13);  // r^2 = 185 *
     350            ADD_RAND ( 74,  8, 11);  // r^2 = 185 *
     351            ADD_RAND ( 75,  7, 12);  // r^2 = 193
     352            ADD_RAND ( 76,  5, 13);  // r^2 = 194
     353            ADD_AXIS ( 77, 14)    ;  // r^2 = 196
     354            if (Ns <= 14) goto finish;
     355
     356            ADD_RAND ( 78,  1, 14);  // r^2 = 197
     357            ADD_DIAG ( 79, 10)    ;  // r^2 = 200 *
     358            ADD_RAND ( 79,  2, 14);  // r^2 = 200 *
     359            ADD_RAND ( 80,  9, 11);  // r^2 = 202
     360            ADD_RAND ( 81,  3, 14);  // r^2 = 205 *
     361            ADD_RAND ( 81,  6, 13);  // r^2 = 205 *
     362            ADD_RAND ( 82,  8, 12);  // r^2 = 208
     363            ADD_RAND ( 83,  4, 14);  // r^2 = 212
     364            ADD_RAND ( 84,  7, 13);  // r^2 = 218
     365            ADD_RAND ( 85, 10, 11);  // r^2 = 221 *
     366            ADD_RAND ( 85,  5, 14);  // r^2 = 221 *
     367            ADD_RAND ( 86,  9, 12);  // r^2 = 225 *
     368            ADD_AXIS ( 86, 15);      // r^2 = 225 *
     369            if (Ns <= 15) goto finish;
     370
    277371        finish:
    278372            vo[iy][ix] = s;
    … …  
    298392}
    299393
     394/*
     395101, 104, 106, 109, 113,
     396116, 117, 121, 122, 125,
     397128, 130, 136, 137, 144,
     398145, 146, 148, 149, 153,
     399157, 160, 162, 164, 169,
     400170, 173, 178, 180, 181,
     401185, 193, 194, 196, 197,
     402200, 202, 205, 208, 212,
     403218, 221, 225
     404*/
     405

  • branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c

    r36282 r36319  
    1414#include "pstap.h"
    1515
    16 bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float sigma, float kappa);
     16bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa);
    1717
    1818int main(int argc, char **argv) {
    … …  
    2323
    2424    // generate an image to convolve
    25     tap_convolve_delta (301, 301, 150, 150, 2.0, 1.0);
    26     tap_convolve_delta (301, 301, 100, 200, 2.0, 1.0);
    27     tap_convolve_delta (301, 301, 200, 100, 2.0, 1.0);
    28     tap_convolve_delta (301, 301, 200, 200, 2.0, 1.0);
    29     tap_convolve_delta (301, 301, 100, 100, 2.0, 1.0);
     25    tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 1.0);
     26    tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 1.0);
     27    tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 1.0);
     28    tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 1.0);
     29    tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 1.0);
    3030
    31     tap_convolve_delta (301, 301, 150, 150, 3.0, 1.0);
    32     tap_convolve_delta (301, 301, 100, 200, 3.0, 1.0);
    33     tap_convolve_delta (301, 301, 200, 100, 3.0, 1.0);
    34     tap_convolve_delta (301, 301, 200, 200, 3.0, 1.0);
    35     tap_convolve_delta (301, 301, 100, 100, 3.0, 1.0);
     31    tap_convolve_delta (301, 301, 150, 150, 3.0, 3.0, 1.0);
     32    tap_convolve_delta (301, 301, 100, 200, 3.0, 3.0, 1.0);
     33    tap_convolve_delta (301, 301, 200, 100, 3.0, 3.0, 1.0);
     34    tap_convolve_delta (301, 301, 200, 200, 3.0, 3.0, 1.0);
     35    tap_convolve_delta (301, 301, 100, 100, 3.0, 3.0, 1.0);
    3636
    37     tap_convolve_delta (301, 301, 150, 150, 1.0, 1.0);
    38     tap_convolve_delta (301, 301, 100, 200, 1.0, 1.0);
    39     tap_convolve_delta (301, 301, 200, 100, 1.0, 1.0);
    40     tap_convolve_delta (301, 301, 200, 200, 1.0, 1.0);
    41     tap_convolve_delta (301, 301, 100, 100, 1.0, 1.0);
     37    tap_convolve_delta (301, 301, 150, 150, 3.0, 1.0, 1.0);
     38    tap_convolve_delta (301, 301, 100, 200, 3.0, 1.0, 1.0);
     39    tap_convolve_delta (301, 301, 200, 100, 3.0, 1.0, 1.0);
     40    tap_convolve_delta (301, 301, 200, 200, 3.0, 1.0, 1.0);
     41    tap_convolve_delta (301, 301, 100, 100, 3.0, 1.0, 1.0);
    4242
    43     tap_convolve_delta (301, 301, 150, 150, 2.0, 2.0);
    44     tap_convolve_delta (301, 301, 100, 200, 2.0, 2.0);
    45     tap_convolve_delta (301, 301, 200, 100, 2.0, 2.0);
    46     tap_convolve_delta (301, 301, 200, 200, 2.0, 2.0);
    47     tap_convolve_delta (301, 301, 100, 100, 2.0, 2.0);
     43    tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 2.0);
     44    tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 2.0);
     45    tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 2.0);
     46    tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 2.0);
     47    tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 2.0);
     48
     49    tap_convolve_delta (301, 301, 150, 150, 4.0, 3.0, 1.0);
     50    tap_convolve_delta (301, 301, 100, 200, 4.0, 3.0, 1.0);
     51    tap_convolve_delta (301, 301, 200, 100, 4.0, 3.0, 1.0);
     52    tap_convolve_delta (301, 301, 200, 200, 4.0, 3.0, 1.0);
     53    tap_convolve_delta (301, 301, 100, 100, 4.0, 3.0, 1.0);
     54
     55    tap_convolve_delta (301, 301, 150, 150, 5.0, 3.0, 1.0);
     56    tap_convolve_delta (301, 301, 100, 200, 5.0, 3.0, 1.0);
     57    tap_convolve_delta (301, 301, 200, 100, 5.0, 3.0, 1.0);
     58    tap_convolve_delta (301, 301, 200, 200, 5.0, 3.0, 1.0);
     59    tap_convolve_delta (301, 301, 100, 100, 5.0, 3.0, 1.0);
     60
     61    tap_convolve_delta (301, 301, 150, 150, 5.0, 1.0, 1.0);
     62    tap_convolve_delta (301, 301, 100, 200, 5.0, 1.0, 1.0);
     63    tap_convolve_delta (301, 301, 200, 100, 5.0, 1.0, 1.0);
     64    tap_convolve_delta (301, 301, 200, 200, 5.0, 1.0, 1.0);
     65    tap_convolve_delta (301, 301, 100, 100, 5.0, 1.0, 1.0);
    4866
    4967    ok(!psMemCheckLeaks (id, NULL, stdout, false), "no memory leaks");
    … …  
    5270}
    5371
    54 bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float sigma, float kappa) {
     72bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa) {
    5573
    5674    bool status;
    … …  
    6078    input->data.F32[Yo][Xo] = 1.0;
    6179
    62     psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (3.0);
     80    psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (Nsigma);
    6381    ok(smdata, "init smdata");
    6482    ok(!smdata->radflux, "init radflux to NULL");
    65     ok(smdata->Nsigma == 3, "init Nsigma to 3");
     83    ok(smdata->Nsigma == Nsigma, "init Nsigma");
    6684   
    6785    status = psImageSmooth2dCacheKernel_PS1_V1 (smdata, sigma, kappa);
    … …  
    7189    ok(status, "smoothed image");
    7290   
    73     float s = 0.0;
     91    double s = 0.0;
    7492
    7593    psImage *model = psImageAlloc (Nx, Ny, PS_TYPE_F32);
    … …  
    7795    for (int iy = 0; iy < model->numRows; iy++) {
    7896        for (int ix = 0; ix < model->numCols; ix++) {
    79             float z = PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + PS_SQR((iy + 0.5 - Yo - 0.5) / sigma);
    80             if (z > PS_SQR(smdata->Nsigma)) continue;
     97            float z = 0.5*PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + 0.5*PS_SQR((iy + 0.5 - Yo - 0.5) / sigma);
     98            if (z > 0.5*PS_SQR(smdata->Nsigma)) continue;
    8199           
    82100            float f = 1.0 / (1.0 + z*kappa + pow(z, 1.666));
    … …  
    101119        }
    102120    }
    103     ok (fabs(min_df) < 1e-5, "minor deviations");
    104     ok (fabs(max_df) < 1e-5, "minor deviations");
     121    ok (fabs(min_df) < 3e-5, "minor deviations");
     122    ok (fabs(max_df) < 3e-5, "minor deviations");
    105123
    106     if ((fabs(min_df) >= 1e-5) || (fabs(max_df) >= 1e-5)) {
     124    if ((fabs(min_df) >= 3e-5) || (fabs(max_df) >= 3e-5)) {
    107125        psFits *fits = NULL;
    108126        fits = psFitsOpen ("input.fits", "w");
Note: See TracChangeset for help on using the changeset viewer.