IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 24, 2005, 3:22:12 PM (21 years ago)
Author:
gusciora
Message:

...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/pmObjects.c

    r3078 r3086  
    55 *  @author GLG, MHPCC
    66 *
    7  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-01-24 22:57:52 $
     7 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-01-25 01:22:12 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2222above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
    2323the location (x value) of all peaks.
     24 
     25XXX: What types should be supported?  Only F32 us implemented.
     26 
     27XXX: We currently step through the input vector twice; once to determine the
     28size of the output vector, then to set the values of the output vector.
     29Depending upon actual use, this may need to be optimized.
    2430 *****************************************************************************/
    2531psVector *pmFindVectorPeeks(const psVector *vector,
    2632                            psF32 threshold)
    2733{
    28     return(NULL);
     34    PS_VECTOR_CHECK_NULL(vector, NULL);
     35    PS_VECTOR_CHECK_TYPE(vector, PS_TYPE_F32, NULL);
     36    int count = 0;
     37    int n = vector->n;
     38
     39    if ((vector->data.F32[0] > vector->data.F32[1]) &&
     40            (vector->data.F32[0] > threshold)) {
     41        count++;
     42    }
     43
     44    for (psU32 i = 1; i < n-1 ; i++) {
     45        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
     46                (vector->data.F32[i] > vector->data.F32[i+1]) &&
     47                (vector->data.F32[i] > threshold)) {
     48            count++;
     49        }
     50    }
     51    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
     52            (vector->data.F32[n-1] > threshold)) {
     53        count++;
     54    }
     55
     56    psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);
     57    count = 0;
     58    if ((vector->data.F32[0] > vector->data.F32[1]) &&
     59            (vector->data.F32[0] > threshold)) {
     60        tmpVector->data.U32[count++] = 0;
     61    }
     62    for (psU32 i = 1; i < (n-1) ; i++) {
     63        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
     64                (vector->data.F32[i] > vector->data.F32[i+1]) &&
     65                (vector->data.F32[i] > threshold)) {
     66            tmpVector->data.U32[count++] = i;
     67        }
     68    }
     69    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
     70            (vector->data.F32[n-1] > threshold)) {
     71        tmpVector->data.U32[count++] = n-1;
     72    }
     73
     74    return(tmpVector);
     75}
     76
     77/******************************************************************************
     78XXX: Is there a better way to do this?
     79 *****************************************************************************/
     80psVector *p_psGetRowVectorFromImage(psImage *image,
     81                                    psU32 row)
     82{
     83    PS_IMAGE_CHECK_NULL(image, NULL);
     84    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
     85
     86    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
     87    for (psU32 col = 0; col < image->numCols ; col++) {
     88        tmpVector->data.F32[col] = image->data.F32[row][col];
     89    }
     90    return(tmpVector);
    2991}
    3092
    3193/******************************************************************************
    3294pmFindImagePeeks(image, threshold): Find all local peaks in the given psImage
    33 above the given threshold.  Returns a psList containing the location (x/y
    34 value) of all peaks.
     95above the given threshold.  Returns a psList containing location (x/y value)
     96of all peaks.
    3597 *****************************************************************************/
    3698psList *pmFindImagePeeks(const psImage *image,
    3799                         psF32 threshold)
    38100{
     101    PS_IMAGE_CHECK_NULL(image, NULL);
     102    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
     103    psVector *tmpRow = NULL;
     104    psU32 col = 0;
     105    psU32 row = 0;
     106    psVector *rowX = NULL;
     107
     108    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 0);
     109    psVector *row1 = pmFindVectorPeeks(tmpRow, threshold);
     110    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 1);
     111    psVector *row2 = pmFindVectorPeeks(tmpRow, threshold);
     112    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 2);
     113    psVector *row3 = pmFindVectorPeeks(tmpRow, threshold);
     114
     115    row = 0;
     116    for (psU32 i = 0 ; i < row1->n ; i++ ) {
     117        col = row1->data.U32[i];
     118        if (col == 0) {
     119            if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     120                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     121                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     122                // Add peak at location (row, col)
     123            }
     124        } else if (col < (image->numCols - 1)) {
     125            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     126                    (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     127                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
     128                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     129                    (image->data.F32[row][col] >= image->data.F32[row+1][col+2])) {
     130                // Add peak at location (row, col)
     131            }
     132
     133        } else if (col == (image->numCols - 1)) {
     134            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     135                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     136                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
     137                // Add peak at location (row, col)
     138            }
     139
     140        } else {
     141            printf("XXX: ERROR!\n");
     142        }
     143    }
     144
     145    row = image->numRows - 1;
     146    for (psU32 i = 0 ; i < row1->n ; i++ ) {
     147        col = row1->data.U32[i];
     148        if (col == 0) {
     149            if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     150                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     151                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
     152                // Add peak at location (row, col)
     153            }
     154        } else if (col < (image->numCols - 1)) {
     155            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     156                    (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
     157                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
     158                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     159                    (image->data.F32[row][col] >= image->data.F32[row+1][col+2])) {
     160                // Add peak at location (row, col)
     161            }
     162
     163        } else if (col == (image->numCols - 1)) {
     164            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
     165                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
     166                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
     167                // Add peak at location (row, col)
     168            }
     169
     170        } else {
     171            printf("XXX: ERROR!\n");
     172        }
     173    }
     174
     175
     176
    39177    return(NULL);
    40178}
    … …  
    141279/******************************************************************************
    142280XXX: Why only *x argument?
     281 
     282May x[0] is x and x[1] is y?
     283 
     284    params->data.F32[0] = So;
     285    params->data.F32[1] = Zo;
     286    params->data.F32[2] = Xo;
     287    params->data.F32[3] = Yo;
     288    params->data.F32[4] = sqrt(2.0) / SigmaX;
     289    params->data.F32[5] = sqrt(2.0) / SigmaY;
     290    params->data.F32[6] = Sxy;
    143291 *****************************************************************************/
    144292psF32 pmMinLM_Gauss2D(psVector *deriv,
    … …  
    146294                      psVector *x)
    147295{
    148     return(0.0);
    149 }
    150 
    151 /******************************************************************************
    152 XXX: Why only *x argument?
     296    psF32 X = x->data.F32[0] - params->data.F32[2];
     297    psF32 Y = x->data.F32[1] - params->data.F32[3];
     298    psF32 px = params->data.F32[4]*X;
     299    psF32 py = params->data.F32[5]*Y;
     300    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     301    psF32 r = exp(-z);
     302    psF32 f = params->data.F32[1]*r + params->data.F32[0];
     303
     304    psF32 q = params->data.F32[1]*r;
     305    deriv->data.F32[0] = +1.0;
     306    deriv->data.F32[1] = +r;
     307    deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
     308    deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
     309    deriv->data.F32[4] = -2.0*q*px*X;
     310    deriv->data.F32[5] = -2.0*q*py*Y;
     311    deriv->data.F32[6] = -q*X*Y;
     312
     313    return(f);
     314}
     315
     316/******************************************************************************
     317    params->data.F32[0] = So;
     318    params->data.F32[1] = Zo;
     319    params->data.F32[2] = Xo;
     320    params->data.F32[3] = Yo;
     321    params->data.F32[4] = sqrt(2) / SigmaX;
     322    params->data.F32[5] = sqrt(2) / SigmaY;
     323    params->data.F32[6] = Sxy;
    153324 *****************************************************************************/
    154325psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
    … …  
    156327                            psVector *x)
    157328{
    158     return(0.0);
    159 }
    160 
    161 /******************************************************************************
     329    psF32 X = x->data.F32[0] - params->data.F32[2];
     330    psF32 Y = x->data.F32[1] - params->data.F32[3];
     331    psF32 px = params->data.F32[4]*X;
     332    psF32 py = params->data.F32[5]*Y;
     333    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     334    psF32 t = 1 + z + 0.5*z*z;
     335    psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
     336    psF32 f = params->data.F32[1]*r + params->data.F32[0];
     337
     338    //
     339    // note difference from a pure gaussian: q = params->data.F32[1]*r
     340    //
     341
     342    psF32 q = params->data.F32[1]*r*r*t;
     343    deriv->data.F32[0] = +1.0;
     344    deriv->data.F32[1] = +r;
     345    deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
     346    deriv->data.F32[3] = q *
     347                         (2.0*py*params->data.F32[5] + params->data.F32[6]*X);
     348    deriv->data.F32[4] = -2.0*q*px*X;
     349    deriv->data.F32[5] = -2.0*q*py*Y;
     350    deriv->data.F32[6] = -q*X*Y;
     351
     352    return(f);
     353}
     354
     355/******************************************************************************
     356    params->data.F32[0] = So;
     357    params->data.F32[1] = Zo;
     358    params->data.F32[2] = Xo;
     359    params->data.F32[3] = Yo;
     360    params->data.F32[4] = Sx;
     361    params->data.F32[5] = Sy;
     362    params->data.F32[6] = Sxy;
     363    params->data.F32[7] = B2;
     364    params->data.F32[8] = B3;
    162365 *****************************************************************************/
    163366psF32 pmMinLM_Wauss2D(psVector *deriv,
    … …  
    165368                      psVector *x)
    166369{
    167     return(0.0);
    168 }
    169 
    170 /******************************************************************************
     370    psF32 X = x->data.F32[0] - params->data.F32[2];
     371    psF32 Y = x->data.F32[1] - params->data.F32[2];
     372    psF32 px = params->data.F32[4]*X;
     373    psF32 py = params->data.F32[5]*Y;
     374    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
     375    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
     376    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
     377    psF32 f = params->data.F32[1]*r + params->data.F32[0];
     378
     379    //
     380    // note difference from gaussian: q = params->data.F32[1]*r
     381    //
     382    psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
     383    deriv->data.F32[0] = +1.0;
     384    deriv->data.F32[1] = +r;
     385    deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
     386    deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
     387    deriv->data.F32[4] = -2.0*q*px*X;
     388    deriv->data.F32[5] = -2.0*q*py*Y;
     389    deriv->data.F32[6] = -q*X*Y;
     390    deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
     391    deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
     392    //
     393    // The values of 100 dampen the swing of params->data.F32[7,8] */
     394    //
     395
     396    return(f);
     397}
     398
     399// XXX: What should these be?
     400#define FFACTOR 1.0
     401#define FSCALE 1.0
     402/******************************************************************************
     403    params->data.F32[0] = So;
     404    params->data.F32[1] = Zo;
     405    params->data.F32[2] = Xo;
     406    params->data.F32[3] = Yo;
     407    params->data.F32[4] = SxInner;
     408    params->data.F32[5] = SyInner;
     409    params->data.F32[6] = SxyInner;
     410    params->data.F32[7] = SxOuter;
     411    params->data.F32[8] = SyOuter;
     412    params->data.F32[9] = SxyOuter;
     413    params->data.F32[10] = N;
    171414 *****************************************************************************/
    172415psF32 pmMinLM_TwistGauss2D(psVector *deriv,
    … …  
    174417                           psVector *x)
    175418{
    176     return(0.0);
    177 }
    178 
    179 /******************************************************************************
     419    psF32 X = x->data.F32[0] - params->data.F32[2];
     420    psF32 Y = x->data.F32[1] - params->data.F32[3];
     421    psF32 px1 = params->data.F32[4]*X;
     422    psF32 py1 = params->data.F32[5]*Y;
     423    psF32 px2 = params->data.F32[7]*X;
     424    psF32 py2 = params->data.F32[8]*Y;
     425    psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
     426    psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
     427    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
     428
     429
     430    psF32 f = params->data.F32[5]*r + params->data.F32[6];
     431    psF32 q1 = params->data.F32[5]*PS_SQR(r);
     432    psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
     433    deriv->data.F32[0] = +1.0;
     434    deriv->data.F32[1] = +r;
     435    deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y);
     436    deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X);
     437
     438    //
     439    // These fudge factors impede the growth of params->data.F32[4] beyond
     440    // params->data.F32[7].
     441    //
     442    psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
     443    psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
     444    deriv->data.F32[4] = -2.0*q1*px1*X*f2;
     445
     446    //
     447    // These fudge factors impede the growth of params->data.F32[5] beyond
     448    // params->data.F32[8].
     449    //
     450    f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
     451    f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
     452    deriv->data.F32[5] = -2.0*q1*py1*Y*f2;
     453    deriv->data.F32[6] = -q1*X*Y;
     454    deriv->data.F32[7] = -2.0*q2*px2*X;
     455    deriv->data.F32[8] = -2.0*q2*py2*Y;
     456    deriv->data.F32[9] = -q2*X*Y;
     457    deriv->data.F32[10] = -q1*log(z2);
     458
     459    return(f);
     460}
     461
     462/******************************************************************************
     463    float Sersic()
     464    params->data.F32[0] = So;
     465    params->data.F32[1] = Zo;
     466    params->data.F32[2] = Xo;
     467    params->data.F32[3] = Yo;
     468    params->data.F32[4] = Sx;
     469    params->data.F32[5] = Sy;
     470    params->data.F32[6] = Sxy;
     471    params->data.F32[7] = Nexp;
    180472 *****************************************************************************/
    181473psF32 pmMinLM_Sersic(psVector *deriv,
    … …  
    187479
    188480/******************************************************************************
     481    float SersicBulge()
     482    params->data.F32[0] So;
     483    params->data.F32[1] Zo;
     484    params->data.F32[2] Xo;
     485    params->data.F32[3] Yo;
     486    params->data.F32[4] SxInner;
     487    params->data.F32[5] SyInner;
     488    params->data.F32[6] SxyInner;
     489    params->data.F32[7] Zd;
     490    params->data.F32[8] SxOuter;
     491    params->data.F32[9] SyOuter;
     492    params->data.F32[10] = SxyOuter;
     493    params->data.F32[11] = Nexp;
    189494 *****************************************************************************/
    190495psF32 pmMinLM_SersicCore(psVector *deriv,
Note: See TracChangeset for help on using the changeset viewer.