IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36375 for trunk/psLib


Ignore:
Timestamp:
Dec 10, 2013, 2:55:11 PM (13 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-20130904

Location:
trunk
Files:
9 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psLib/src/imageops/Makefile.am

    r21207 r36375  
    77        psImageBackground.c \
    88        psImageConvolve.c \
     9        psImageConvolve2dCache.c \
    910        psImageCovariance.c \
    1011        psImageGeomManip.c \

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

    r35767 r36375  
    3636#define MIN_GAUSS_FRAC 0.25             // Minimum Gaussian fraction to accept when smoothing
    3737
    38 
    39 
    4038static bool threaded = false;           // Run image convolution threaded?
    4139static pthread_mutex_t threadMutex = PTHREAD_MUTEX_INITIALIZER;
    42 
    4340
    4441static void kernelFree(psKernel *kernel)
    … …  
    780777}
    781778
    782 void psImageSmooth_PreAlloc_DataFree (psImageSmooth_PreAlloc_Data *smdata) {
     779void psImageSmoothCacheDataFree (psImageSmoothCacheData *smdata) {
    783780    psFree (smdata->resultX);
    784781    psFree (smdata->resultY);
    … …  
    786783}
    787784
    788 psImageSmooth_PreAlloc_Data *psImageSmooth_PreAlloc_DataAlloc (psImage *image, double sigma, double Nsigma) {
    789 
    790     psImageSmooth_PreAlloc_Data *smdata = psAlloc(sizeof(psImageSmooth_PreAlloc_Data));
    791     psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth_PreAlloc_DataFree);
     785// allocate the psImageSmoothCache data structure, but do not define the kernel
     786psImageSmoothCacheData *psImageSmoothCacheAlloc (psImage *image, double sigma, double Nsigma) {
     787
     788    psImageSmoothCacheData *smdata = psAlloc(sizeof(psImageSmoothCacheData));
     789    psMemSetDeallocator(smdata, (psFreeFunc) psImageSmoothCacheDataFree);
     790
     791    smdata->kernel = NULL;
    792792
    793793    if (!image) {
    … …  
    796796        smdata->Nx = 0;
    797797        smdata->Ny = 0;
    798         smdata->kernel = NULL;
    799798        smdata->resultX = NULL;
    800799        smdata->resultY = NULL;
    … …  
    807806    smdata->Ny = image->numRows;            // Number of rows
    808807
    809     IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32);
     808    // XXX drop this : we now require a call to a kernel-creation function (like psImageSmoothCacheKernel_Gauss)
     809    // IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32);
    810810       
    811811    // use a temp running buffer for X and Y directions.
    … …  
    819819}
    820820
     821// generate a Gaussian smoothing kernel for supplied sigma.  sigma here does not need to match
     822// that used to allocate the structure, but it is recommended
     823bool psImageSmoothCacheKernel_Gauss (psImageSmoothCacheData *smdata, float sigma) {
     824    // check for NULL structure elements?
     825    psFree (smdata->kernel);
     826    IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32);
     827    return true;
     828}
     829
    821830// we can use the same DATA structure on multiple images of the same size
    822 bool psImageSmooth_PreAlloc_F32(psImage *image, psImageSmooth_PreAlloc_Data *smdata)
     831bool psImageSmoothCache_F32(psImage *image, psImageSmoothCacheData *smdata)
    823832{
    824833    PS_ASSERT_IMAGE_NON_NULL(image, false);
     834    PS_ASSERT_VECTOR_NON_NULL(smdata->kernel, false);
    825835    // assert on data type
    826836
    … …  
    21562166    return threaded;
    21572167}
    2158 

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

    r35767 r36375  
    2626#define PS_TYPE_KERNEL_NAME "psF32"    ///< the data type for kernel as a string */
    2727
    28 /// a structure to contain data related to image smoothing with a 1D gauss kernel
     28/// a structure to contain data related to image smoothing with a pre-cached 1D gauss kernel
    2929typedef struct {
    3030    int Nx;
    … …  
    3434    psF32 *resultY;
    3535    psVector *kernel;
    36 } psImageSmooth_PreAlloc_Data;
     36} psImageSmoothCacheData;
     37
     38/// a structure to contain data related to image smoothing with a pre-cached 1D gauss kernel
     39typedef struct {
     40    float Nsigma;
     41    int Ns;                             // number of pixel radii
     42    float *radflux;                     // conv kernel in special positions
     43} psImageSmooth2dCacheData;
    3744
    3845/// A convolution kernel
    … …  
    306313);
    307314
    308 psImageSmooth_PreAlloc_Data *psImageSmooth_PreAlloc_DataAlloc (psImage *image, double sigma, double Nsigma);
    309 bool psImageSmooth_PreAlloc_F32(psImage *image, psImageSmooth_PreAlloc_Data *smdata);
     315psImageSmoothCacheData *psImageSmoothCacheAlloc (psImage *image, double sigma, double Nsigma);
     316bool psImageSmoothCache_F32(psImage *image, psImageSmoothCacheData *smdata);
     317bool psImageSmoothCacheKernel_Gauss (psImageSmoothCacheData *smdata, float sigma);
    310318
    311319/// Control threading for image convolution functions
    … …  
    318326bool psImageConvolveGetThreads(void);
    319327
     328psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma);
     329bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa);
     330bool psImageSmooth2dCacheKernel_Gauss (psImageSmooth2dCacheData *smdata, float sigma);
     331bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata);
     332
    320333/// @}
    321334#endif // #ifndef PS_IMAGE_CONVOLVE_H

  • trunk/psLib/src/math/psEllipse.c

    r35555 r36375  
    1212// f = exp(-z) where z describes the elliptical contour at any flux level:
    1313
     14// NOTE: major, minor are the 1-sigma lengths in the major,minor directions assuming the
     15// moments represent a Gaussian profile
     16
    1417// sigma shape: z = 0.5((x/sx)^2 + (y/sy)^2 + sxy*x*y)
    1518// sigma axes : z = 0.5((x/sa)^2 + (y/sb)^2), x,y rotated by theta
    … …  
    1720// polarization : e0, e1, e2
    1821
    19 // ellipse rotation (major, minor, theta) -> (x2, y2, xy)
     22// ellipse rotation (major, minor, theta) -> (Mxx, Mxy, Myy)
     23psEllipseMoments psEllipseAxesToMoments(psEllipseAxes axes)
     24{
     25    psEllipseMoments moments;
     26
     27    double f1 = PS_SQR(axes.major) + PS_SQR(axes.minor);
     28    double f2 = PS_SQR(axes.major) - PS_SQR(axes.minor);
     29
     30    moments.x2 = +0.5*f1 + 0.5*f2*cos(2*axes.theta);
     31    moments.y2 = +0.5*f1 - 0.5*f2*cos(2*axes.theta);
     32    moments.xy = +0.5*f2*sin(2*axes.theta);
     33
     34    assert (isfinite(moments.x2));
     35    assert (isfinite(moments.y2));
     36    assert (isfinite(moments.xy));
     37
     38    return moments;
     39}
     40
     41// ellipse rotation (Mxx, Mxy, Myy) -> (major, minor, theta).
     42psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments, double maxAR)
     43{
     44    psEllipseAxes axes;
     45    psEllipseAxes badValue = {NAN, NAN, NAN};
     46
     47    if (!isfinite(moments.x2)) return badValue;
     48    if (!isfinite(moments.y2)) return badValue;
     49    if (!isfinite(moments.xy)) return badValue;
     50
     51    if (moments.x2 < 0) return badValue;
     52    if (moments.y2 < 0) return badValue;
     53
     54    double g1 = moments.x2 + moments.y2;
     55    double g2 = moments.x2 - moments.y2;
     56    double g3 = sqrt(PS_SQR(g2) + 4*PS_SQR(moments.xy));
     57
     58    axes.major = sqrt (0.5*(g1 + g3));
     59    axes.theta = +0.5 * atan2 (+2.0*moments.xy, g2); // theta in radians
     60
     61    // long, thin objects are likely to have a poorly measured minor axis
     62    // the angle and major axis are likely to be ok.
     63    // restrict the axis ratio
     64    double rAR2 = (g1 - g3) / (g1 + g3);
     65    if (rAR2 < 1.0/PS_SQR(maxAR)) {
     66        axes.minor = axes.major / maxAR;
     67    } else {
     68        axes.minor = sqrt (0.5*(g1 - g3));
     69    }
     70
     71    assert (isfinite(axes.major));
     72    assert (isfinite(axes.minor));
     73    assert (isfinite(axes.theta));
     74
     75    return axes;
     76}
     77
     78// ellipse rotation (major, minor, theta) -> (sx, sy, sxy)
     79// theta is postive rotation of major axis away from x-axis
     80psEllipseShape psEllipseAxesToShape(psEllipseAxes axes)
     81{
     82    psEllipseShape shape;
     83    psEllipseShape badValue = {NAN, NAN, NAN};
     84
     85    if (!isfinite(axes.minor)) return badValue;
     86    if (!isfinite(axes.major)) return badValue;
     87    if (!isfinite(axes.theta)) return badValue;
     88
     89    if (axes.minor <= 0) return badValue;
     90    if (axes.major <= 0) return badValue;
     91
     92    double f1 = 1.0 / PS_SQR(axes.minor) + 1.0 / PS_SQR(axes.major);
     93    double f2 = 1.0 / PS_SQR(axes.minor) - 1.0 / PS_SQR(axes.major);
     94
     95    double sxr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
     96    double syr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
     97
     98    // sxr, syr cannot be < 0 (f1 >= f2)
     99
     100    shape.sx  = +1.0 / sqrt(sxr);
     101    shape.sy  = +1.0 / sqrt(syr);
     102    shape.sxy = -0.5*f2*sin(2*axes.theta);
     103
     104    assert (isfinite(shape.sx));
     105    assert (isfinite(shape.sy));
     106    assert (isfinite(shape.sxy));
     107
     108    return shape;
     109}
     110
     111// ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
     112psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape, double maxAR)
     113{
     114    psEllipseAxes axes;
     115
     116    double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx);
     117    double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx);
     118    double f3 = sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
     119
     120    axes.minor = sqrt (2.0 / (f1 + f3));
     121    axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2);
     122
     123    // long, thin objects are likely to have a poorly measured major axis
     124    // the angle and minor axis are likely to be ok.
     125    // restrict the axis ratio
     126    double rAR2 = (f1 - f3) / (f1 + f3);
     127    if (rAR2 < 1.0/PS_SQR(maxAR)) {
     128        axes.major = axes.minor * maxAR;
     129    } else {
     130        axes.major = sqrt (2.0 / (f1 - f3));
     131    }
     132
     133    assert (isfinite(axes.theta));
     134    assert (isfinite(axes.major));
     135    assert (isfinite(axes.minor));
     136
     137    return axes;
     138}
     139
     140// ellipse rotation (major, minor, theta) -> (e0, e1, e2)
    20141psEllipsePol psEllipseAxesToPol(psEllipseAxes axes)
    21142{
    … …  
    35156}
    36157
    37 // ellipse rotation (major, minor, theta) -> (x2, y2, xy)
    38 psEllipsePol psEllipseShapeToPol(psEllipseShape shape)
    39 {
    40     psEllipsePol pol;
    41 
    42     double r = 1.0 / (1.0 - PS_SQR(shape.sx)*PS_SQR(shape.sy)*PS_SQR(shape.sxy));
    43 
    44     pol.e0 = r*(PS_SQR(shape.sx) + PS_SQR(shape.sy));
    45     pol.e1 = r*(PS_SQR(shape.sx) - PS_SQR(shape.sy));
    46     // XXX I do not understand this negative sign
    47     pol.e2 = -r*(2.0*PS_SQR(shape.sx)*PS_SQR(shape.sy)*shape.sxy);
    48 
    49     assert (isfinite(pol.e0));
    50     assert (isfinite(pol.e1));
    51     assert (isfinite(pol.e2));
    52 
    53     return pol;
    54 }
    55 
    56 // ellipse rotation (major, minor, theta) -> (x2, y2, xy)
    57 // XXXX handle case where e0 < LIMIT
     158// ellipse rotation (e0, e1, e2) -> (major, minor, theta)
    58159psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol,
    59160                                 const float minMinorAxis)
    … …  
    96197}
    97198
    98 // ellipse rotation (major, minor, theta) -> (x2, y2, xy)
    99 psEllipseMoments psEllipseAxesToMoments(psEllipseAxes axes)
    100 {
    101     psEllipseMoments moments;
    102 
    103     double f1 = PS_SQR(axes.major) + PS_SQR(axes.minor);
    104     double f2 = PS_SQR(axes.major) - PS_SQR(axes.minor);
    105 
    106     moments.x2 = +0.5*f1 + 0.5*f2*cos(2*axes.theta);
    107     moments.y2 = +0.5*f1 - 0.5*f2*cos(2*axes.theta);
    108     moments.xy = +0.5*f2*sin(2*axes.theta);
    109 
    110     assert (isfinite(moments.x2));
    111     assert (isfinite(moments.y2));
    112     assert (isfinite(moments.xy));
    113 
    114     return moments;
    115 }
    116 
    117 // ellipse rotation (x2, y2, xy) -> (major, minor, theta).  NOTE: major, minor are the
    118 // 1-sigma lengths in the major,minor directions assuming the moments represent a Gaussian
    119 // profile
    120 psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments, double maxAR)
    121 {
    122     psEllipseAxes axes;
    123     psEllipseAxes badValue = {NAN, NAN, NAN};
    124 
    125     if (!isfinite(moments.x2)) return badValue;
    126     if (!isfinite(moments.y2)) return badValue;
    127     if (!isfinite(moments.xy)) return badValue;
    128 
    129     if (moments.x2 < 0) return badValue;
    130     if (moments.y2 < 0) return badValue;
    131 
    132     double g1 = moments.x2 + moments.y2;
    133     double g2 = moments.x2 - moments.y2;
    134     double g3 = sqrt(PS_SQR(g2) + 4*PS_SQR(moments.xy));
    135 
    136     axes.major = sqrt (0.5*(g1 + g3));
    137     axes.theta = +0.5 * atan2 (+2.0*moments.xy, g2); // theta in radians
    138 
    139     // long, thin objects are likely to have a poorly measured minor axis
    140     // the angle and major axis are likely to be ok.
    141     // restrict the axis ratio
    142     double rAR2 = (g1 - g3) / (g1 + g3);
    143     if (rAR2 < 1.0/PS_SQR(maxAR)) {
    144         axes.minor = axes.major / maxAR;
    145     } else {
    146         axes.minor = sqrt (0.5*(g1 - g3));
    147     }
    148 
    149     assert (isfinite(axes.major));
    150     assert (isfinite(axes.minor));
    151     assert (isfinite(axes.theta));
    152 
    153     return axes;
    154 }
    155 
    156 // ellipse rotation (major, minor, theta) -> (sx, sy, sxy)
    157 // theta is postive rotation of major axis away from x-axis
    158 psEllipseShape psEllipseAxesToShape(psEllipseAxes axes)
     199// ellipse rotation (sx, sy, sxy) -> (e0, e1, e2)
     200psEllipsePol psEllipseShapeToPol(psEllipseShape shape)
     201{
     202    psEllipsePol pol;
     203
     204    double r = 1.0 / (1.0 - PS_SQR(shape.sx)*PS_SQR(shape.sy)*PS_SQR(shape.sxy));
     205
     206    pol.e0 = r*(PS_SQR(shape.sx) + PS_SQR(shape.sy));
     207    pol.e1 = r*(PS_SQR(shape.sx) - PS_SQR(shape.sy));
     208    // XXX I do not understand this negative sign
     209    pol.e2 = -r*(2.0*PS_SQR(shape.sx)*PS_SQR(shape.sy)*shape.sxy);
     210
     211    assert (isfinite(pol.e0));
     212    assert (isfinite(pol.e1));
     213    assert (isfinite(pol.e2));
     214
     215    return pol;
     216}
     217
     218// ellipse rotation (e0, e1, e2) -> (sx, sy, sxy)
     219psEllipseShape psEllipsePolToShape(psEllipsePol pol)
    159220{
    160221    psEllipseShape shape;
    161     psEllipseShape badValue = {NAN, NAN, NAN};
    162 
    163     if (!isfinite(axes.minor)) return badValue;
    164     if (!isfinite(axes.major)) return badValue;
    165     if (!isfinite(axes.theta)) return badValue;
    166 
    167     if (axes.minor <= 0) return badValue;
    168     if (axes.major <= 0) return badValue;
    169 
    170     double f1 = 1.0 / PS_SQR(axes.minor) + 1.0 / PS_SQR(axes.major);
    171     double f2 = 1.0 / PS_SQR(axes.minor) - 1.0 / PS_SQR(axes.major);
    172 
    173     double sxr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
    174     double syr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
    175 
    176     // sxr, syr cannot be < 0 (f1 >= f2)
    177 
    178     shape.sx  = +1.0 / sqrt(sxr);
    179     shape.sy  = +1.0 / sqrt(syr);
    180     shape.sxy = -0.5*f2*sin(2*axes.theta);
     222
     223    double q = sqrt (PS_SQR(pol.e1) + PS_SQR(pol.e2));
     224    double p = PS_SQR(pol.e0) + PS_SQR(q) - 2.0*q*pol.e0;
     225
     226    // double f1 = 4*pol.e0 / p;
     227    // double f2 = 4*q      / p;
     228   
     229    double sxr = 2.0*(pol.e0 - pol.e1) / p;
     230    double syr = 2.0*(pol.e0 + pol.e1) / p;
     231
     232    shape.sx = sqrt(1.0 / sxr);
     233    shape.sy = sqrt(1.0 / syr);
     234
     235    shape.sxy = -2.0 * pol.e2 / p;
    181236
    182237    assert (isfinite(shape.sx));
    … …  
    187242}
    188243
    189 // ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
    190 psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape, double maxAR)
    191 {
    192     psEllipseAxes axes;
    193 
    194     double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx);
    195     double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx);
    196     double f3 = sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
    197 
    198     axes.minor = sqrt (2.0 / (f1 + f3));
    199     axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2);
    200 
    201     // long, thin objects are likely to have a poorly measured major axis
    202     // the angle and minor axis are likely to be ok.
    203     // restrict the axis ratio
    204     double rAR2 = (f1 - f3) / (f1 + f3);
    205     if (rAR2 < 1.0/PS_SQR(maxAR)) {
    206         axes.major = axes.minor * maxAR;
    207     } else {
    208         axes.major = sqrt (2.0 / (f1 - f3));
    209     }
    210 
    211     assert (isfinite(axes.theta));
    212     assert (isfinite(axes.major));
    213     assert (isfinite(axes.minor));
    214 
    215     return axes;
    216 }
    217 
    218244// XXX keep this construction?
    219245// force the axis ratio to be less than 10

  • trunk/psLib/src/math/psEllipse.h

    r19996 r36375  
    3838} psEllipseShape;
    3939
    40 /// Ellipse defined in terms of polarisations
     40/// Ellipse defined in terms of polarizations
    4141typedef struct {
    4242    double e0;                          ///< Scale (Mxx + Myy)
    … …  
    6565                                   );
    6666
    67 /// Convert axes to polarisation representation
     67/// Convert axes to polarization representation
    6868psEllipsePol psEllipseAxesToPol(psEllipseAxes axes ///< Axes of ellipse
    6969                                );
    7070
    71 /// Convert shape to polarisation representation
     71/// Convert polarization to axes representation
     72psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol, ///< Polarization of ellipse
     73                                 const float minMinorAxis ///< Minimum allowed minor axis
     74    );
     75
     76/// Convert shape to polarization representation
    7277psEllipsePol psEllipseShapeToPol(psEllipseShape shape ///< Shape of ellipse
    7378                                 );
    7479
    75 /// Convert polarisation to axes representation
    76 ///
    77 /// XXX This API goes against the PS convention of outputs being first.
    78 psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol, ///< Polarisation of ellipse
    79                                  const float minMinorAxis ///< Minimum allowed minor axis
    80     );
     80/// Convert shape to polarization representation
     81psEllipseShape psEllipsePolToShape(psEllipsePol pol ///< Shape of ellipse
     82                                 );
    8183
    8284/// @}

  • trunk/psLib/src/math/psMixtureModels.c

  • trunk/psLib/test/imageops

    • Property svn:ignore

      •  

        old new  
        2525tap_psImageConvolve
        2626tap_psImageConvolve2
         27tap_psImageConvolve2dCache
        2728convolutionBench
        2829tap_psImageInterpolate2

  • trunk/psLib/test/imageops/Makefile.am

    r35767 r36375  
    2222        tap_psImageConvolve \
    2323        tap_psImageConvolve2 \
     24        tap_psImageConvolve2dCache \
    2425        tap_psImagePixelExtract \
    2526        tap_psImageInterpolate2 \
Note: See TracChangeset for help on using the changeset viewer.