IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 22, 2006, 11:12:20 AM (20 years ago)
Author:
magnier
Message:

fixed minor bugs in pmSourceFitModelInit

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_rel9_p0/psModules/src/objects/pmObjects.c

    r6182 r6183  
    66 *  @author EAM, IfA: significant modifications.
    77 *
    8  *  @version $Revision: 1.5.4.3 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-01-22 20:20:47 $
     8 *  @version $Revision: 1.5.4.4 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-01-22 21:12:19 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    884884            // XXX EAM : should this limit be user-defined?
    885885            if (PS_SQR(pDiff) < wDiff) {
    886 
    887                 Var += wDiff;
    888                 Sum += pDiff;
    889 
    890                 psF32 xWght = xDiff * pDiff;
    891                 psF32 yWght = yDiff * pDiff;
    892 
    893                 X1  += xWght;
    894                 Y1  += yWght;
    895                 XY  += xDiff * yWght;
    896 
    897                 X2  += xDiff * xWght;
    898                 Y2  += yDiff * yWght;
    899 
    900                 peakPixel = PS_MAX (vPix[row][col], peakPixel);
    901                 numPixels++;
    902             }
    903         }
    904 
    905         // if we have less than (1/4) of the possible pixels, force a retry
    906         // XXX EAM - the limit is a bit arbitrary.  make it user defined?
    907         if ((numPixels < 0.75*R2) || (Sum <= 0)) {
    908             psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
    909             psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    910             return (false);
    911         }
    912 
    913         psTrace (".psModules.pmSourceMoments", 5,
    914                  "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
    915                  sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
    916 
    917         //
    918         // first moment X  = X1/Sum + xc
    919         // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
    920         // Sxy             = XY / Sum
    921         //
    922         x = X1/Sum;
    923         y = Y1/Sum;
    924         if ((fabs(x) > radius) || (fabs(y) > radius)) {
    925             psTrace (".psModules.pmSourceMoments", 3,
    926                      "large centroid swing; invalid peak %d, %d\n",
    927                      source->peak->x, source->peak->y);
    928             psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    929             return (false);
    930         }
    931 
    932         source->moments->x = x + xPeak;
    933         source->moments->y = y + yPeak;
    934 
    935         // XXX EAM : Sxy needs to have x*y subtracted
    936         source->moments->Sxy = XY/Sum - x*y;
    937         source->moments->Sum = Sum;
    938         source->moments->SN  = Sum / sqrt(Var);
    939         source->moments->Peak = peakPixel;
    940         source->moments->nPixels = numPixels;
    941 
    942         // XXX EAM : these values can be negative, so we need to limit the range
    943         source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
    944         source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
    945 
    946         psTrace (".psModules.pmSourceMoments", 4,
    947                  "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
    948                  sky, Sum, source->moments->x, source->moments->y,
    949                  source->moments->Sx, source->moments->Sy, source->moments->Sxy);
    950 
    951         psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    952         return(true);
    953     }
    954 
    955     // XXX EAM : I used
    956     int pmComparePeakAscend (const void **a, const void **b) {
    957         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    958         pmPeak *A = *(pmPeak **)a;
    959         pmPeak *B = *(pmPeak **)b;
    960 
    961         psF32 diff;
    962 
    963         diff = A->counts - B->counts;
    964         if (diff < FLT_EPSILON) {
    965             psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
    966             return (-1);
    967         } else if (diff > FLT_EPSILON) {
    968             psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
    969             return (+1);
    970         }
    971         psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
    972         return (0);
    973     }
    974 
    975     int pmComparePeakDescend (const void **a, const void **b) {
    976         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    977         pmPeak *A = *(pmPeak **)a;
    978         pmPeak *B = *(pmPeak **)b;
    979 
    980         psF32 diff;
    981 
    982         diff = A->counts - B->counts;
    983         if (diff < FLT_EPSILON) {
    984             psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
    985             return (+1);
    986         } else if (diff > FLT_EPSILON) {
    987             psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
    988             return (-1);
    989         }
    990         psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
    991         return (0);
    992     }
    993 
    994     /******************************************************************************
     886                continue;
     887            }
     888
     889            Var += wDiff;
     890            Sum += pDiff;
     891
     892            psF32 xWght = xDiff * pDiff;
     893            psF32 yWght = yDiff * pDiff;
     894
     895            X1  += xWght;
     896            Y1  += yWght;
     897            XY  += xDiff * yWght;
     898
     899            X2  += xDiff * xWght;
     900            Y2  += yDiff * yWght;
     901
     902            peakPixel = PS_MAX (vPix[row][col], peakPixel);
     903            numPixels++;
     904        }
     905    }
     906
     907    // if we have less than (1/4) of the possible pixels, force a retry
     908    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     909    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
     910        psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n");
     911        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     912        return (false);
     913    }
     914
     915    psTrace (".psModules.pmSourceMoments", 5,
     916             "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n",
     917             sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
     918
     919    //
     920    // first moment X  = X1/Sum + xc
     921    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
     922    // Sxy             = XY / Sum
     923    //
     924    x = X1/Sum;
     925    y = Y1/Sum;
     926    if ((fabs(x) > radius) || (fabs(y) > radius)) {
     927        psTrace (".psModules.pmSourceMoments", 3,
     928                 "large centroid swing; invalid peak %d, %d\n",
     929                 source->peak->x, source->peak->y);
     930        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     931        return (false);
     932    }
     933
     934    source->moments->x = x + xPeak;
     935    source->moments->y = y + yPeak;
     936
     937    // XXX EAM : Sxy needs to have x*y subtracted
     938    source->moments->Sxy = XY/Sum - x*y;
     939    source->moments->Sum = Sum;
     940    source->moments->SN  = Sum / sqrt(Var);
     941    source->moments->Peak = peakPixel;
     942    source->moments->nPixels = numPixels;
     943
     944    // XXX EAM : these values can be negative, so we need to limit the range
     945    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
     946    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
     947
     948    psTrace (".psModules.pmSourceMoments", 4,
     949             "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n",
     950             sky, Sum, source->moments->x, source->moments->y,
     951             source->moments->Sx, source->moments->Sy, source->moments->Sxy);
     952
     953    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     954    return(true);
     955}
     956
     957// XXX EAM : I used
     958int pmComparePeakAscend (const void **a, const void **b)
     959{
     960    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     961    pmPeak *A = *(pmPeak **)a;
     962    pmPeak *B = *(pmPeak **)b;
     963
     964    psF32 diff;
     965
     966    diff = A->counts - B->counts;
     967    if (diff < FLT_EPSILON) {
     968        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
     969        return (-1);
     970    } else if (diff > FLT_EPSILON) {
     971        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
     972        return (+1);
     973    }
     974    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     975    return (0);
     976}
     977
     978int pmComparePeakDescend (const void **a, const void **b)
     979{
     980    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     981    pmPeak *A = *(pmPeak **)a;
     982    pmPeak *B = *(pmPeak **)b;
     983
     984    psF32 diff;
     985
     986    diff = A->counts - B->counts;
     987    if (diff < FLT_EPSILON) {
     988        psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__);
     989        return (+1);
     990    } else if (diff > FLT_EPSILON) {
     991        psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__);
     992        return (-1);
     993    }
     994    psTrace(__func__, 3, "---- %s(0) end ----\n", __func__);
     995    return (0);
     996}
     997
     998/******************************************************************************
    995999    pmSourcePSFClump(source, metadata): Find the likely PSF clump in the
    9961000    sigma-x, sigma-y plane. return 0,0 clump in case of error.
    997     *****************************************************************************/
    998 
    999     // XXX EAM include a S/N cutoff in selecting the sources?
    1000     pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) {
    1001         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1002 
    1003         # define NPIX 10
    1004         # define SCALE 0.1
    1005 
    1006         psArray *peaks  = NULL;
    1007         pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
    1008         pmPSFClump psfClump = emptyClump;
    1009 
    1010         PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
    1011         PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
    1012 
    1013         // find the sigmaX, sigmaY clump
     1001*****************************************************************************/
     1002
     1003// XXX EAM include a S/N cutoff in selecting the sources?
     1004pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
     1005{
     1006    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1007
     1008    # define NPIX 10
     1009    # define SCALE 0.1
     1010
     1011    psArray *peaks  = NULL;
     1012    pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0};
     1013    pmPSFClump psfClump = emptyClump;
     1014
     1015    PS_ASSERT_PTR_NON_NULL(sources, emptyClump);
     1016    PS_ASSERT_PTR_NON_NULL(metadata, emptyClump);
     1017
     1018    // find the sigmaX, sigmaY clump
     1019    {
     1020        psStats *stats  = NULL;
     1021        psImage *splane = NULL;
     1022        int binX, binY;
     1023        bool status;
     1024
     1025        psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
     1026        if (!status)
     1027            SX_MAX = 10.0;
     1028        psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
     1029        if (!status)
     1030            SY_MAX = 10.0;
     1031
     1032        // construct a sigma-plane image
     1033        // psImageAlloc does zero the data
     1034        splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
     1035        for (int i = 0; i < splane->numRows; i++)
    10141036        {
    1015             psStats *stats  = NULL;
    1016             psImage *splane = NULL;
    1017             int binX, binY;
    1018             bool status;
    1019 
    1020             psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX");
    1021             if (!status)
    1022                 SX_MAX = 10.0;
    1023             psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX");
    1024             if (!status)
    1025                 SY_MAX = 10.0;
    1026 
    1027             // construct a sigma-plane image
    1028             // psImageAlloc does zero the data
    1029             splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32);
    1030             for (int i = 0; i < splane->numRows; i++)
    1031             {
    1032                 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
    1033             }
    1034 
    1035             // place the sources in the sigma-plane image (ignore 0,0 values?)
    1036             for (psS32 i = 0 ; i < sources->n ; i++)
    1037             {
    1038                 pmSource *tmpSrc = (pmSource *) sources->data[i];
    1039                 if (tmpSrc == NULL) {
    1040                     continue;
    1041                 }
    1042                 if (tmpSrc->moments == NULL) {
    1043                     continue;
    1044                 }
    1045 
    1046                 // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
    1047                 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
    1048                     continue;
    1049                 }
    1050 
    1051                 // for the moment, force splane dimensions to be 10x10 image pix
    1052                 binX = tmpSrc->moments->Sx/SCALE;
    1053                 if (binX < 0)
    1054                     continue;
    1055                 if (binX >= splane->numCols)
    1056                     continue;
    1057 
    1058                 binY = tmpSrc->moments->Sy/SCALE;
    1059                 if (binY < 0)
    1060                     continue;
    1061                 if (binY >= splane->numRows)
    1062                     continue;
    1063 
    1064                 splane->data.F32[binY][binX] += 1.0;
    1065             }
    1066 
    1067             // find the peak in this image
    1068             stats = psStatsAlloc (PS_STAT_MAX);
    1069             stats = psImageStats (stats, splane, NULL, 0);
    1070             peaks = pmFindImagePeaks (splane, stats[0].max / 2);
    1071             psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
    1072 
    1073             psFree (splane);
    1074             psFree (stats);
    1075 
    1076         }
    1077         // XXX EAM : possible errors:
    1078         //           1) no peak in splane
    1079         //           2) no significant peak in splane
    1080 
    1081         // measure statistics on Sx, Sy if Sx, Sy within range of clump
     1037            memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32));
     1038        }
     1039
     1040        // place the sources in the sigma-plane image (ignore 0,0 values?)
     1041        for (psS32 i = 0 ; i < sources->n ; i++)
    10821042        {
    1083             pmPeak *clump;
    1084             psF32 minSx, maxSx;
    1085             psF32 minSy, maxSy;
    1086             psVector *tmpSx = NULL;
    1087             psVector *tmpSy = NULL;
    1088             psStats *stats  = NULL;
    1089 
    1090             // XXX EAM : this lets us takes the single highest peak
    1091             psArraySort (peaks, pmComparePeakDescend);
    1092             clump = peaks->data[0];
    1093             psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
    1094 
    1095             // define section window for clump
    1096             minSx = clump->x * SCALE - 0.2;
    1097             maxSx = clump->x * SCALE + 0.2;
    1098             minSy = clump->y * SCALE - 0.2;
    1099             maxSy = clump->y * SCALE + 0.2;
    1100 
    1101             tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
    1102             tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
    1103             tmpSx->n = 0;
    1104             tmpSy->n = 0;
    1105 
    1106             // XXX clip sources based on flux?
    1107             // create vectors with Sx, Sy values in window
    1108             for (psS32 i = 0 ; i < sources->n ; i++)
    1109             {
    1110                 pmSource *tmpSrc = (pmSource *) sources->data[i];
    1111 
    1112                 if (tmpSrc->moments->Sx < minSx)
    1113                     continue;
    1114                 if (tmpSrc->moments->Sx > maxSx)
    1115                     continue;
    1116                 if (tmpSrc->moments->Sy < minSy)
    1117                     continue;
    1118                 if (tmpSrc->moments->Sy > maxSy)
    1119                     continue;
    1120                 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
    1121                 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
    1122                 tmpSx->n++;
    1123                 tmpSy->n++;
    1124                 if (tmpSx->n == tmpSx->nalloc) {
    1125                     psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
    1126                     psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
    1127                 }
    1128             }
    1129 
    1130             // measures stats of Sx, Sy
    1131             stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    1132 
    1133             stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
    1134             psfClump.X  = stats->clippedMean;
    1135             psfClump.dX = stats->clippedStdev;
    1136 
    1137             stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
    1138             psfClump.Y  = stats->clippedMean;
    1139             psfClump.dY = stats->clippedStdev;
    1140 
    1141             psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
    1142             psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
    1143             // these values should be pushed on the metadata somewhere
    1144 
    1145             psFree (stats);
    1146             psFree (peaks);
    1147             psFree (tmpSx);
    1148             psFree (tmpSy);
    1149         }
    1150 
    1151         psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    1152         return (psfClump);
    1153     }
    1154 
    1155     /******************************************************************************
     1043            pmSource *tmpSrc = (pmSource *) sources->data[i];
     1044            if (tmpSrc == NULL) {
     1045                continue;
     1046            }
     1047            if (tmpSrc->moments == NULL) {
     1048                continue;
     1049            }
     1050
     1051            // Sx,Sy are limited at 0.  a peak at 0,0 is artificial
     1052            if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) {
     1053                continue;
     1054            }
     1055
     1056            // for the moment, force splane dimensions to be 10x10 image pix
     1057            binX = tmpSrc->moments->Sx/SCALE;
     1058            if (binX < 0)
     1059                continue;
     1060            if (binX >= splane->numCols)
     1061                continue;
     1062
     1063            binY = tmpSrc->moments->Sy/SCALE;
     1064            if (binY < 0)
     1065                continue;
     1066            if (binY >= splane->numRows)
     1067                continue;
     1068
     1069            splane->data.F32[binY][binX] += 1.0;
     1070        }
     1071
     1072        // find the peak in this image
     1073        stats = psStatsAlloc (PS_STAT_MAX);
     1074        stats = psImageStats (stats, splane, NULL, 0);
     1075        peaks = pmFindImagePeaks (splane, stats[0].max / 2);
     1076        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2);
     1077
     1078        psFree (splane);
     1079        psFree (stats);
     1080
     1081    }
     1082    // XXX EAM : possible errors:
     1083    //           1) no peak in splane
     1084    //           2) no significant peak in splane
     1085
     1086    // measure statistics on Sx, Sy if Sx, Sy within range of clump
     1087    {
     1088        pmPeak *clump;
     1089        psF32 minSx, maxSx;
     1090        psF32 minSy, maxSy;
     1091        psVector *tmpSx = NULL;
     1092        psVector *tmpSy = NULL;
     1093        psStats *stats  = NULL;
     1094
     1095        // XXX EAM : this lets us takes the single highest peak
     1096        psArraySort (peaks, pmComparePeakDescend);
     1097        clump = peaks->data[0];
     1098        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts);
     1099
     1100        // define section window for clump
     1101        minSx = clump->x * SCALE - 0.2;
     1102        maxSx = clump->x * SCALE + 0.2;
     1103        minSy = clump->y * SCALE - 0.2;
     1104        maxSy = clump->y * SCALE + 0.2;
     1105
     1106        tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32);
     1107        tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32);
     1108        tmpSx->n = 0;
     1109        tmpSy->n = 0;
     1110
     1111        // XXX clip sources based on flux?
     1112        // create vectors with Sx, Sy values in window
     1113        for (psS32 i = 0 ; i < sources->n ; i++)
     1114        {
     1115            pmSource *tmpSrc = (pmSource *) sources->data[i];
     1116
     1117            if (tmpSrc->moments->Sx < minSx)
     1118                continue;
     1119            if (tmpSrc->moments->Sx > maxSx)
     1120                continue;
     1121            if (tmpSrc->moments->Sy < minSy)
     1122                continue;
     1123            if (tmpSrc->moments->Sy > maxSy)
     1124                continue;
     1125            tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx;
     1126            tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy;
     1127            tmpSx->n++;
     1128            tmpSy->n++;
     1129            if (tmpSx->n == tmpSx->nalloc) {
     1130                psVectorRealloc (tmpSx, tmpSx->nalloc + 100);
     1131                psVectorRealloc (tmpSy, tmpSy->nalloc + 100);
     1132            }
     1133        }
     1134
     1135        // measures stats of Sx, Sy
     1136        stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     1137
     1138        stats = psVectorStats (stats, tmpSx, NULL, NULL, 0);
     1139        psfClump.X  = stats->clippedMean;
     1140        psfClump.dX = stats->clippedStdev;
     1141
     1142        stats = psVectorStats (stats, tmpSy, NULL, NULL, 0);
     1143        psfClump.Y  = stats->clippedMean;
     1144        psfClump.dY = stats->clippedStdev;
     1145
     1146        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump  X,  Y: %f, %f\n", psfClump.X, psfClump.Y);
     1147        psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
     1148        // these values should be pushed on the metadata somewhere
     1149
     1150        psFree (stats);
     1151        psFree (peaks);
     1152        psFree (tmpSx);
     1153        psFree (tmpSy);
     1154    }
     1155
     1156    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     1157    return (psfClump);
     1158}
     1159
     1160/******************************************************************************
    11561161    pmSourceRoughClass(source, metadata): make a guess at the source
    11571162    classification.
    … …  
    11621167     
    11631168    EAM: I moved S/N calculation to pmSourceMoments, using weight image
    1164     *****************************************************************************/
    1165 
    1166     bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) {
    1167         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1168 
    1169         psBool rc = true;
    1170 
    1171         int Nsat     = 0;
    1172         int Next     = 0;
    1173         int Nstar    = 0;
    1174         int Npsf     = 0;
    1175         int Ncr      = 0;
    1176         int Nsatstar = 0;
    1177         // psRegion allArray = psRegionSet (0, 0, 0, 0);
    1178         psRegion inner;
    1179 
    1180         // report stats on S/N values for star-like objects
    1181         psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
    1182         starsn->n = 0;
    1183 
    1184         // check return status value (do these exist?)
    1185         bool status;
    1186         psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
    1187 
    1188         // XXX allow clump size to be scaled relative to sigmas?
    1189         // make rough IDs based on clumpX,Y,DX,DY
    1190         for (psS32 i = 0 ; i < sources->n ; i++) {
    1191 
    1192             pmSource *tmpSrc = (pmSource *) sources->data[i];
    1193 
    1194             tmpSrc->peak->class = 0;
    1195 
    1196             psF32 sigX = tmpSrc->moments->Sx;
    1197             psF32 sigY = tmpSrc->moments->Sy;
    1198 
    1199             // XXX EAM : can we use the value of SATURATE if mask is NULL?
    1200             inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
    1201             int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
    1202 
    1203             // saturated star (size consistent with PSF or larger)
    1204             // Nsigma should be user-configured parameter
    1205             bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
    1206             big = true;
    1207             if ((Nsatpix > 1) && big) {
    1208                 tmpSrc->type = PM_SOURCE_STAR;
    1209                 tmpSrc->mode = PM_SOURCE_SATSTAR;
    1210                 Nsatstar ++;
    1211                 continue;
    1212             }
    1213 
    1214             // saturated object (not a star, eg bleed trails, hot pixels)
    1215             if (Nsatpix > 1) {
    1216                 tmpSrc->type = PM_SOURCE_SATURATED;
    1217                 tmpSrc->mode = PM_SOURCE_DEFAULT;
    1218                 Nsat ++;
    1219                 continue;
    1220             }
    1221 
    1222             // likely defect (too small to be stellar) (push out to 3 sigma)
    1223             // low S/N objects which are small are probably stellar
    1224             // only set candidate defects if
    1225             if ((sigX < 0.05) || (sigY < 0.05)) {
    1226                 tmpSrc->type = PM_SOURCE_DEFECT;
    1227                 tmpSrc->mode = PM_SOURCE_DEFAULT;
    1228                 Ncr ++;
    1229                 continue;
    1230             }
    1231 
    1232             // likely unsaturated extended source (too large to be stellar)
    1233             if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
    1234                 tmpSrc->type = PM_SOURCE_EXTENDED;
    1235                 tmpSrc->mode = PM_SOURCE_DEFAULT;
    1236                 Next ++;
    1237                 continue;
    1238             }
    1239 
    1240             // the rest are probable stellar objects
    1241             starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
    1242             starsn->n ++;
    1243             Nstar ++;
    1244 
    1245             // PSF star (within 1.5 sigma of clump center, S/N > limit)
    1246             psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
    1247             if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
    1248                 tmpSrc->type = PM_SOURCE_STAR;
    1249                 tmpSrc->mode = PM_SOURCE_PSFSTAR;
    1250                 Npsf ++;
    1251                 continue;
    1252             }
    1253 
    1254             // random type of star
     1169*****************************************************************************/
     1170
     1171bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
     1172{
     1173    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1174
     1175    psBool rc = true;
     1176
     1177    int Nsat     = 0;
     1178    int Next     = 0;
     1179    int Nstar    = 0;
     1180    int Npsf     = 0;
     1181    int Ncr      = 0;
     1182    int Nsatstar = 0;
     1183    // psRegion allArray = psRegionSet (0, 0, 0, 0);
     1184    psRegion inner;
     1185
     1186    // report stats on S/N values for star-like objects
     1187    psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32);
     1188    starsn->n = 0;
     1189
     1190    // check return status value (do these exist?)
     1191    bool status;
     1192    psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM");
     1193
     1194    // XXX allow clump size to be scaled relative to sigmas?
     1195    // make rough IDs based on clumpX,Y,DX,DY
     1196    for (psS32 i = 0 ; i < sources->n ; i++) {
     1197
     1198        pmSource *tmpSrc = (pmSource *) sources->data[i];
     1199
     1200        tmpSrc->peak->class = 0;
     1201
     1202        psF32 sigX = tmpSrc->moments->Sx;
     1203        psF32 sigY = tmpSrc->moments->Sy;
     1204
     1205        // XXX EAM : can we use the value of SATURATE if mask is NULL?
     1206        inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2);
     1207        int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED);
     1208
     1209        // saturated star (size consistent with PSF or larger)
     1210        // Nsigma should be user-configured parameter
     1211        bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
     1212        big = true;
     1213        if ((Nsatpix > 1) && big) {
    12551214            tmpSrc->type = PM_SOURCE_STAR;
     1215            tmpSrc->mode = PM_SOURCE_SATSTAR;
     1216            Nsatstar ++;
     1217            continue;
     1218        }
     1219
     1220        // saturated object (not a star, eg bleed trails, hot pixels)
     1221        if (Nsatpix > 1) {
     1222            tmpSrc->type = PM_SOURCE_SATURATED;
    12561223            tmpSrc->mode = PM_SOURCE_DEFAULT;
    1257         }
    1258 
    1259         {
    1260             psStats *stats  = NULL;
    1261             stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
    1262             stats = psVectorStats (stats, starsn, NULL, NULL, 0);
    1263             psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
    1264             psFree (stats);
    1265             psFree (starsn);
    1266         }
    1267 
    1268         psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
    1269         psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
    1270         psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
    1271         psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
    1272         psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
    1273         psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
    1274 
    1275         psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1276         return(rc);
    1277     }
    1278 
    1279     /** pmSourceDefinePixels()
    1280      *
    1281      * Define psImage subarrays for the source located at coordinates x,y on the
    1282      * image set defined by readout. The pixels defined by this operation consist of
    1283      * a square window (of full width 2Radius+1) centered on the pixel which contains
    1284      * the given coordinate, in the frame of the readout. The window is defined to
    1285      * have limits which are valid within the boundary of the readout image, thus if
    1286      * the radius would fall outside the image pixels, the subimage is truncated to
    1287      * only consist of valid pixels. If readout->mask or readout->weight are not
    1288      * NULL, matching subimages are defined for those images as well. This function
    1289      * fails if no valid pixels can be defined (x or y less than Radius, for
    1290      * example). This function should be used to define a region of interest around a
    1291      * source, including both source and sky pixels.
    1292      *
    1293      * XXX: must code this.
    1294      *
    1295      */
    1296     bool pmSourceDefinePixels(
    1297         pmSource *mySource,                 ///< Add comment.
    1298         pmReadout *readout,                 ///< Add comment.
    1299         psF32 x,                            ///< Add comment.
    1300         psF32 y,                            ///< Add comment.
    1301         psF32 Radius)                       ///< Add comment.
     1224            Nsat ++;
     1225            continue;
     1226        }
     1227
     1228        // likely defect (too small to be stellar) (push out to 3 sigma)
     1229        // low S/N objects which are small are probably stellar
     1230        // only set candidate defects if
     1231        if ((sigX < 0.05) || (sigY < 0.05)) {
     1232            tmpSrc->type = PM_SOURCE_DEFECT;
     1233            tmpSrc->mode = PM_SOURCE_DEFAULT;
     1234            Ncr ++;
     1235            continue;
     1236        }
     1237
     1238        // likely unsaturated extended source (too large to be stellar)
     1239        if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
     1240            tmpSrc->type = PM_SOURCE_EXTENDED;
     1241            tmpSrc->mode = PM_SOURCE_DEFAULT;
     1242            Next ++;
     1243            continue;
     1244        }
     1245
     1246        // the rest are probable stellar objects
     1247        starsn->data.F32[starsn->n] = tmpSrc->moments->SN;
     1248        starsn->n ++;
     1249        Nstar ++;
     1250
     1251        // PSF star (within 1.5 sigma of clump center, S/N > limit)
     1252        psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
     1253        if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) {
     1254            tmpSrc->type = PM_SOURCE_STAR;
     1255            tmpSrc->mode = PM_SOURCE_PSFSTAR;
     1256            Npsf ++;
     1257            continue;
     1258        }
     1259
     1260        // random type of star
     1261        tmpSrc->type = PM_SOURCE_STAR;
     1262        tmpSrc->mode = PM_SOURCE_DEFAULT;
     1263    }
     1264
    13021265    {
    1303         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1304         psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented.  Returning FALSE.\n");
    1305         psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    1306         return(false);
    1307     }
    1308 
    1309     /******************************************************************************
     1266        psStats *stats  = NULL;
     1267        stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
     1268        stats = psVectorStats (stats, starsn, NULL, NULL, 0);
     1269        psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max);
     1270        psFree (stats);
     1271        psFree (starsn);
     1272    }
     1273
     1274    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar:    %3d\n", Nstar);
     1275    psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf:     %3d\n", Npsf);
     1276    psTrace (".pmObjects.pmSourceRoughClass", 2, "Next:     %3d\n", Next);
     1277    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar);
     1278    psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat:     %3d\n", Nsat);
     1279    psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr:      %3d\n", Ncr);
     1280
     1281    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     1282    return(rc);
     1283}
     1284
     1285/** pmSourceDefinePixels()
     1286 *
     1287 * Define psImage subarrays for the source located at coordinates x,y on the
     1288 * image set defined by readout. The pixels defined by this operation consist of
     1289 * a square window (of full width 2Radius+1) centered on the pixel which contains
     1290 * the given coordinate, in the frame of the readout. The window is defined to
     1291 * have limits which are valid within the boundary of the readout image, thus if
     1292 * the radius would fall outside the image pixels, the subimage is truncated to
     1293 * only consist of valid pixels. If readout->mask or readout->weight are not
     1294 * NULL, matching subimages are defined for those images as well. This function
     1295 * fails if no valid pixels can be defined (x or y less than Radius, for
     1296 * example). This function should be used to define a region of interest around a
     1297 * source, including both source and sky pixels.
     1298 *
     1299 * XXX: must code this.
     1300 *
     1301 */
     1302bool pmSourceDefinePixels(
     1303    pmSource *mySource,                 ///< Add comment.
     1304    pmReadout *readout,                 ///< Add comment.
     1305    psF32 x,                            ///< Add comment.
     1306    psF32 y,                            ///< Add comment.
     1307    psF32 Radius)                       ///< Add comment.
     1308{
     1309    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1310    psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented.  Returning FALSE.\n");
     1311    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     1312    return(false);
     1313}
     1314
     1315/******************************************************************************
    13101316    pmSourceSetPixelsCircle(source, image, radius)
    13111317     
    13121318    XXX: This was replaced by DefinePixels in SDRS.  Remove it.
    1313     *****************************************************************************/
    1314     bool pmSourceSetPixelsCircle(pmSource *source,
    1315                                  const psImage *image,
    1316                                  psF32 radius) {
    1317         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1318         PS_ASSERT_IMAGE_NON_NULL(image, false);
    1319         PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    1320         PS_ASSERT_PTR_NON_NULL(source, false);
    1321         PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1322         PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1323         PS_FLOAT_COMPARE(0.0, radius, false);
    1324 
    1325         //
    1326         // We define variables for code readability.
    1327         //
    1328         // XXX: Since the peak->xy coords are in image, not subImage coords,
    1329         // these variables should be renamed for clarity (imageCenterRow, etc).
    1330         //
    1331         psS32 radiusS32 = (psS32) radius;
    1332         psS32 SubImageCenterRow = source->peak->y;
    1333         psS32 SubImageCenterCol = source->peak->x;
    1334         // XXX EAM : for the circle to stay on the image
    1335         // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
    1336         psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
    1337         psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
    1338         psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
    1339         psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
    1340 
    1341         // XXX: Must recycle image.
    1342         // XXX EAM: this message reflects a programming error we know about.
    1343         //          i am setting it to a trace message which we can take out
    1344         if (source->pixels != NULL) {
    1345             psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,
    1346                      "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
    1347             psFree(source->pixels);
    1348         }
    1349         source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
    1350                                        SubImageStartRow,
    1351                                        SubImageEndCol,
    1352                                        SubImageEndRow));
    1353 
    1354         // XXX: Must recycle image.
    1355         if (source->mask != NULL) {
    1356             psFree(source->mask);
    1357         }
    1358         source->mask = psImageAlloc(source->pixels->numCols,
    1359                                     source->pixels->numRows,
    1360                                     PS_TYPE_U8); // XXX EAM : type was F32
    1361 
    1362         //
    1363         // Loop through the subimage mask, initialize mask to 0 or 1.
    1364         // XXX EAM: valid pixels should have 0, not 1
    1365         for (psS32 row = 0 ; row < source->mask->numRows; row++) {
    1366             for (psS32 col = 0 ; col < source->mask->numCols; col++) {
    1367 
    1368                 if (checkRadius2((psF32) radiusS32,
    1369                                  (psF32) radiusS32,
    1370                                  radius,
    1371                                  (psF32) col,
    1372                                  (psF32) row)) {
    1373                     source->mask->data.U8[row][col] = 0;
    1374                 } else {
    1375                     source->mask->data.U8[row][col] = 1;
    1376                 }
    1377             }
    1378         }
    1379         psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    1380         return(true);
    1381     }
    1382 
    1383     /******************************************************************************
     1319*****************************************************************************/
     1320bool pmSourceSetPixelsCircle(pmSource *source,
     1321                             const psImage *image,
     1322                             psF32 radius)
     1323{
     1324    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1325    PS_ASSERT_IMAGE_NON_NULL(image, false);
     1326    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
     1327    PS_ASSERT_PTR_NON_NULL(source, false);
     1328    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1329    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1330    PS_FLOAT_COMPARE(0.0, radius, false);
     1331
     1332    //
     1333    // We define variables for code readability.
     1334    //
     1335    // XXX: Since the peak->xy coords are in image, not subImage coords,
     1336    // these variables should be renamed for clarity (imageCenterRow, etc).
     1337    //
     1338    psS32 radiusS32 = (psS32) radius;
     1339    psS32 SubImageCenterRow = source->peak->y;
     1340    psS32 SubImageCenterCol = source->peak->x;
     1341    // XXX EAM : for the circle to stay on the image
     1342    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
     1343    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - radiusS32);
     1344    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1);
     1345    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - radiusS32);
     1346    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1);
     1347
     1348    // XXX: Must recycle image.
     1349    // XXX EAM: this message reflects a programming error we know about.
     1350    //          i am setting it to a trace message which we can take out
     1351    if (source->pixels != NULL) {
     1352        psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4,
     1353                 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL.  Freeing and reallocating.\n");
     1354        psFree(source->pixels);
     1355    }
     1356    source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
     1357                                   SubImageStartRow,
     1358                                   SubImageEndCol,
     1359                                   SubImageEndRow));
     1360
     1361    // XXX: Must recycle image.
     1362    if (source->mask != NULL) {
     1363        psFree(source->mask);
     1364    }
     1365    source->mask = psImageAlloc(source->pixels->numCols,
     1366                                source->pixels->numRows,
     1367                                PS_TYPE_U8); // XXX EAM : type was F32
     1368
     1369    //
     1370    // Loop through the subimage mask, initialize mask to 0 or 1.
     1371    // XXX EAM: valid pixels should have 0, not 1
     1372    for (psS32 row = 0 ; row < source->mask->numRows; row++) {
     1373        for (psS32 col = 0 ; col < source->mask->numCols; col++) {
     1374
     1375            if (checkRadius2((psF32) radiusS32,
     1376                             (psF32) radiusS32,
     1377                             radius,
     1378                             (psF32) col,
     1379                             (psF32) row)) {
     1380                source->mask->data.U8[row][col] = 0;
     1381            } else {
     1382                source->mask->data.U8[row][col] = 1;
     1383            }
     1384        }
     1385    }
     1386    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     1387    return(true);
     1388}
     1389
     1390/******************************************************************************
    13841391    pmSourceModelGuess(source, model): This function allocates a new
    13851392    pmModel structure based on the given modelType specified in the argument list. 
    … …  
    13901397    image, not subImage coords.  Therefore, the calls to the model evaluation
    13911398    functions will be in image, not subImage coords.  Remember this.
    1392     *****************************************************************************/
    1393     pmModel *pmSourceModelGuess(pmSource *source,
    1394                                 pmModelType modelType) {
    1395         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1396         PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1397         PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1398 
    1399         pmModel *model = pmModelAlloc(modelType);
    1400 
    1401         pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
    1402         modelGuessFunc(model, source);
    1403         psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    1404         return(model);
    1405     }
    1406 
    1407     /******************************************************************************
     1399*****************************************************************************/
     1400pmModel *pmSourceModelGuess(pmSource *source,
     1401                            pmModelType modelType)
     1402{
     1403    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1404    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1405    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1406
     1407    pmModel *model = pmModelAlloc(modelType);
     1408
     1409    pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType);
     1410    modelGuessFunc(model, source);
     1411    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     1412    return(model);
     1413}
     1414
     1415/******************************************************************************
    14081416    evalModel(source, level, row): a private function which evaluates the
    14091417    source->modelPSF function at the specified coords.  The coords are subImage, not
    … …  
    14251433    XXX: For a while, the first psVectorAlloc() was generating a seg fault during
    14261434    testing.  Try to reproduce that and debug.
    1427     *****************************************************************************/
    1428 
    1429     // XXX EAM : I have made this a public function
    1430     // XXX EAM : this now uses a pmModel as the input
    1431     // XXX EAM : it was using src->type to find the model, not model->type
    1432     psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) {
    1433         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1434         PS_ASSERT_PTR_NON_NULL(image, false);
    1435         PS_ASSERT_PTR_NON_NULL(model, false);
    1436         PS_ASSERT_PTR_NON_NULL(model->params, false);
    1437 
    1438         // Allocate the x coordinate structure and convert row/col to image space.
    1439         //
    1440         psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1441         x->data.F32[0] = (psF32) (col + image->col0);
    1442         x->data.F32[1] = (psF32) (row + image->row0);
    1443         psF32 tmpF;
    1444         pmModelFunc modelFunc;
    1445 
    1446         modelFunc = pmModelFunc_GetFunction (model->type);
    1447         tmpF = modelFunc (NULL, model->params, x);
    1448         psFree(x);
    1449         psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    1450         return(tmpF);
    1451     }
    1452 
    1453     /******************************************************************************
     1435*****************************************************************************/
     1436
     1437// XXX EAM : I have made this a public function
     1438// XXX EAM : this now uses a pmModel as the input
     1439// XXX EAM : it was using src->type to find the model, not model->type
     1440psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row)
     1441{
     1442    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1443    PS_ASSERT_PTR_NON_NULL(image, false);
     1444    PS_ASSERT_PTR_NON_NULL(model, false);
     1445    PS_ASSERT_PTR_NON_NULL(model->params, false);
     1446
     1447    // Allocate the x coordinate structure and convert row/col to image space.
     1448    //
     1449    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     1450    x->data.F32[0] = (psF32) (col + image->col0);
     1451    x->data.F32[1] = (psF32) (row + image->row0);
     1452    psF32 tmpF;
     1453    pmModelFunc modelFunc;
     1454
     1455    modelFunc = pmModelFunc_GetFunction (model->type);
     1456    tmpF = modelFunc (NULL, model->params, x);
     1457    psFree(x);
     1458    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     1459    return(tmpF);
     1460}
     1461
     1462/******************************************************************************
    14541463    pmSourceContour(src, img, level, mode): For an input subImage, and model, this
    14551464    routine returns a psArray of coordinates that evaluate to the specified level.
    … …  
    14631472    XXX EAM : this function is using the model for the contour, but it should
    14641473              be using only the image counts
    1465     *****************************************************************************/
    1466     psArray *pmSourceContour(pmSource *source,
    1467                              const psImage *image,
    1468                              psF32 level,
    1469                              pmContourType mode) {
    1470         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1471         PS_ASSERT_PTR_NON_NULL(source, false);
    1472         PS_ASSERT_PTR_NON_NULL(image, false);
    1473         PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1474         PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1475         PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1476         PS_ASSERT_PTR_NON_NULL(source->modelEXT, false);
    1477         // XXX EAM : what is the purpose of modelPSF/modelEXT?
    1478 
    1479         //
    1480         // Allocate data for x/y pairs.
    1481         //
    1482         psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
    1483         psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
    1484 
    1485         //
    1486         // Start at the row with peak pixel, then decrement.
    1487         //
    1488         psS32 col = source->peak->x;
    1489         for (psS32 row = source->peak->y; row>= 0 ; row--) {
    1490             // XXX: yVec contain no real information.  Do we really need it?
    1491             yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
    1492             yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
    1493 
    1494             // Starting at peak pixel, search leftwards for the column intercept.
    1495             psF32 leftIntercept = findValue(source, level, row, col, 0);
    1496             if (isnan(leftIntercept)) {
    1497                 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    1498                 psFree(xVec);
    1499                 psFree(yVec);
    1500                 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
    1501                 return(NULL);
    1502                 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
    1503             }
    1504             xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
    1505 
    1506             // Starting at peak pixel, search rightwards for the column intercept.
    1507 
    1508             psF32 rightIntercept = findValue(source, level, row, col, 1);
    1509             if (isnan(rightIntercept)) {
    1510                 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    1511                 psFree(xVec);
    1512                 psFree(yVec);
    1513                 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
    1514                 return(NULL);
    1515                 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
    1516             }
    1517             psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
    1518             xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
    1519 
    1520             // Set starting column for next row
    1521             col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
    1522         }
    1523         //
    1524         // Start at the row (+1) with peak pixel, then increment.
    1525         //
    1526         col = source->peak->x;
    1527         for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {
    1528             // XXX: yVec contain no real information.  Do we really need it?
    1529             yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
    1530             yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
    1531 
    1532             // Starting at peak pixel, search leftwards for the column intercept.
    1533             psF32 leftIntercept = findValue(source, level, row, col, 0);
    1534             if (isnan(leftIntercept)) {
    1535                 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    1536                 psFree(xVec);
    1537                 psFree(yVec);
    1538                 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
    1539                 return(NULL);
    1540                 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
    1541             }
    1542             xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
    1543 
    1544             // Starting at peak pixel, search rightwards for the column intercept.
    1545             psF32 rightIntercept = findValue(source, level, row, col, 1);
    1546             if (isnan(rightIntercept)) {
    1547                 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
    1548                 psFree(xVec);
    1549                 psFree(yVec);
    1550                 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
    1551                 return(NULL);
    1552                 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
    1553             }
    1554             xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
    1555 
    1556             // Set starting column for next row
    1557             col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
    1558         }
    1559 
    1560         //
    1561         // Allocate an array for result, store coord vectors there.
    1562         //
    1563         psArray *tmpArray = psArrayAlloc(2);
    1564         tmpArray->data[0] = (psPtr *) yVec;
    1565         tmpArray->data[1] = (psPtr *) xVec;
    1566         psTrace(__func__, 3, "---- %s() end ----\n", __func__);
    1567         return(tmpArray);
    1568     }
    1569 
    1570     // save a static values so they may be set externally
    1571     static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
    1572     static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
    1573 
    1574     bool pmSourceFitModelInit (float nIter, float tol) {
    1575 
    1576         PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
    1577         PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
    1578     }
    1579 
    1580     bool pmSourceFitModel (pmSource *source,
    1581                            pmModel *model,
    1582                            const bool PSF) {
    1583         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1584         PS_ASSERT_PTR_NON_NULL(source, false);
    1585         PS_ASSERT_PTR_NON_NULL(source->moments, false);
    1586         PS_ASSERT_PTR_NON_NULL(source->peak, false);
    1587         PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    1588         PS_ASSERT_PTR_NON_NULL(source->mask, false);
    1589         PS_ASSERT_PTR_NON_NULL(source->weight, false);
    1590 
    1591         // XXX EAM : is it necessary for the mask & weight to exist?  the
    1592         //           tests below could be conditions (!NULL)
    1593 
    1594         psBool fitStatus = true;
    1595         psBool onPic     = true;
    1596         psBool rc        = true;
    1597 
    1598         psVector *params = model->params;
    1599         psVector *dparams = model->dparams;
    1600         psVector *paramMask = NULL;
    1601 
    1602         pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1603 
    1604         int nParams = PSF ? 4 : params->n;
    1605 
    1606         // maximum number of valid pixels
    1607         psS32 nPix = source->pixels->numRows * source->pixels->numCols;
    1608 
    1609         // construct the coordinate and value entries
    1610         psArray *x = psArrayAlloc(nPix);
    1611         psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
    1612         psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
    1613 
    1614         nPix = 0;
    1615         for (psS32 i = 0; i < source->pixels->numRows; i++) {
    1616             for (psS32 j = 0; j < source->pixels->numCols; j++) {
    1617                 if (source->mask->data.U8[i][j]) {
    1618                     continue;
    1619                 }
    1620                 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    1621 
    1622                 // Convert i/j to image space:
    1623                 coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    1624                 coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    1625                 x->data[nPix] = (psPtr *) coord;
    1626                 y->data.F32[nPix] = source->pixels->data.F32[i][j];
    1627 
    1628                 // psMinimizeLMChi2 takes wt = 1/dY^2
    1629                 if (source->weight->data.F32[i][j] == 0) {
    1630                     continue;
    1631                 }
    1632                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
    1633                 nPix++;
    1634             }
    1635         }
    1636         if (nPix <  nParams + 1) {
    1637             psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
    1638             psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
    1639             model->status = PM_MODEL_BADARGS;
    1640             psFree (x);
    1641             psFree (y);
    1642             psFree (yErr);
    1643             return(false);
    1644         }
    1645         x->n = nPix;
    1646         y->n = nPix;
    1647         yErr->n = nPix;
    1648 
    1649         psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    1650                                 PM_SOURCE_FIT_MODEL_TOLERANCE);
    1651 
    1652         // PSF model only fits first 4 parameters, EXT model fits all
    1653         if (PSF) {
    1654             paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    1655             for (int i = 0; i < 4; i++) {
    1656                 paramMask->data.U8[i] = 0;
    1657             }
    1658             for (int i = 4; i < paramMask->n; i++) {
    1659                 paramMask->data.U8[i] = 1;
    1660             }
    1661         }
    1662 
    1663         // Set the parameter range checks
    1664         pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    1665         psVector *beta_lim = NULL;
    1666         psVector *params_min = NULL;
    1667         psVector *params_max = NULL;
    1668 
    1669         // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
    1670         //           in the SDRS, I define a new psMinimization which will take these in
    1671         psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
    1672         modelLimits (&beta_lim, &params_min, &params_max);
    1673         for (int i = 0; i < params->n; i++) {
    1674             covar->data.F64[0][i] = beta_lim->data.F32[i];
    1675             covar->data.F64[1][i] = params_min->data.F32[i];
    1676             covar->data.F64[2][i] = params_max->data.F32[i];
    1677         }
    1678 
    1679         psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
    1680         fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     1474*****************************************************************************/
     1475psArray *pmSourceContour(pmSource *source,
     1476                         const psImage *image,
     1477                         psF32 level,
     1478                         pmContourType mode)
     1479{
     1480    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1481    PS_ASSERT_PTR_NON_NULL(source, false);
     1482    PS_ASSERT_PTR_NON_NULL(image, false);
     1483    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1484    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1485    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     1486    PS_ASSERT_PTR_NON_NULL(source->modelEXT, false);
     1487    // XXX EAM : what is the purpose of modelPSF/modelEXT?
     1488
     1489    //
     1490    // Allocate data for x/y pairs.
     1491    //
     1492    psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
     1493    psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32);
     1494
     1495    //
     1496    // Start at the row with peak pixel, then decrement.
     1497    //
     1498    psS32 col = source->peak->x;
     1499    for (psS32 row = source->peak->y; row>= 0 ; row--) {
     1500        // XXX: yVec contain no real information.  Do we really need it?
     1501        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
     1502        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
     1503
     1504        // Starting at peak pixel, search leftwards for the column intercept.
     1505        psF32 leftIntercept = findValue(source, level, row, col, 0);
     1506        if (isnan(leftIntercept)) {
     1507            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
     1508            psFree(xVec);
     1509            psFree(yVec);
     1510            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
     1511            return(NULL);
     1512            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
     1513        }
     1514        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
     1515
     1516        // Starting at peak pixel, search rightwards for the column intercept.
     1517
     1518        psF32 rightIntercept = findValue(source, level, row, col, 1);
     1519        if (isnan(rightIntercept)) {
     1520            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
     1521            psFree(xVec);
     1522            psFree(yVec);
     1523            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
     1524            return(NULL);
     1525            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
     1526        }
     1527        psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept);
     1528        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
     1529
     1530        // Set starting column for next row
     1531        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
     1532    }
     1533    //
     1534    // Start at the row (+1) with peak pixel, then increment.
     1535    //
     1536    col = source->peak->x;
     1537    for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) {
     1538        // XXX: yVec contain no real information.  Do we really need it?
     1539        yVec->data.F32[row] = (psF32) (source->pixels->row0 + row);
     1540        yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row);
     1541
     1542        // Starting at peak pixel, search leftwards for the column intercept.
     1543        psF32 leftIntercept = findValue(source, level, row, col, 0);
     1544        if (isnan(leftIntercept)) {
     1545            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
     1546            psFree(xVec);
     1547            psFree(yVec);
     1548            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
     1549            return(NULL);
     1550            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
     1551        }
     1552        xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept;
     1553
     1554        // Starting at peak pixel, search rightwards for the column intercept.
     1555        psF32 rightIntercept = findValue(source, level, row, col, 1);
     1556        if (isnan(rightIntercept)) {
     1557            psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)");
     1558            psFree(xVec);
     1559            psFree(yVec);
     1560            psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__);
     1561            return(NULL);
     1562            //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n");
     1563        }
     1564        xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept;
     1565
     1566        // Set starting column for next row
     1567        col = (psS32) ((leftIntercept + rightIntercept) / 2.0);
     1568    }
     1569
     1570    //
     1571    // Allocate an array for result, store coord vectors there.
     1572    //
     1573    psArray *tmpArray = psArrayAlloc(2);
     1574    tmpArray->data[0] = (psPtr *) yVec;
     1575    tmpArray->data[1] = (psPtr *) xVec;
     1576    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
     1577    return(tmpArray);
     1578}
     1579
     1580// save a static values so they may be set externally
     1581static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
     1582static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
     1583
     1584bool pmSourceFitModelInit (float nIter, float tol)
     1585{
     1586
     1587    PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter;
     1588    PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
     1589    return true;
     1590}
     1591
     1592bool pmSourceFitModel (pmSource *source,
     1593                       pmModel *model,
     1594                       const bool PSF)
     1595{
     1596    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1597    PS_ASSERT_PTR_NON_NULL(source, false);
     1598    PS_ASSERT_PTR_NON_NULL(source->moments, false);
     1599    PS_ASSERT_PTR_NON_NULL(source->peak, false);
     1600    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
     1601    PS_ASSERT_PTR_NON_NULL(source->mask, false);
     1602    PS_ASSERT_PTR_NON_NULL(source->weight, false);
     1603
     1604    // XXX EAM : is it necessary for the mask & weight to exist?  the
     1605    //           tests below could be conditions (!NULL)
     1606
     1607    psBool fitStatus = true;
     1608    psBool onPic     = true;
     1609    psBool rc        = true;
     1610
     1611    psVector *params = model->params;
     1612    psVector *dparams = model->dparams;
     1613    psVector *paramMask = NULL;
     1614
     1615    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     1616
     1617    int nParams = PSF ? 4 : params->n;
     1618
     1619    // maximum number of valid pixels
     1620    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     1621
     1622    // construct the coordinate and value entries
     1623    psArray *x = psArrayAlloc(nPix);
     1624    psVector *y = psVectorAlloc(nPix, PS_TYPE_F32);
     1625    psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32);
     1626
     1627    nPix = 0;
     1628    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     1629        for (psS32 j = 0; j < source->pixels->numCols; j++) {
     1630            if (source->mask->data.U8[i][j]) {
     1631                continue;
     1632            }
     1633            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     1634
     1635            // Convert i/j to image space:
     1636            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     1637            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     1638            x->data[nPix] = (psPtr *) coord;
     1639            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     1640
     1641            // psMinimizeLMChi2 takes wt = 1/dY^2
     1642            if (source->weight->data.F32[i][j] == 0) {
     1643                continue;
     1644            }
     1645            yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     1646            nPix++;
     1647        }
     1648    }
     1649    if (nPix <  nParams + 1) {
     1650        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
     1651        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     1652        model->status = PM_MODEL_BADARGS;
     1653        psFree (x);
     1654        psFree (y);
     1655        psFree (yErr);
     1656        return(false);
     1657    }
     1658    x->n = nPix;
     1659    y->n = nPix;
     1660    yErr->n = nPix;
     1661
     1662    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
     1663                            PM_SOURCE_FIT_MODEL_TOLERANCE);
     1664
     1665    // PSF model only fits first 4 parameters, EXT model fits all
     1666    if (PSF) {
     1667        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     1668        for (int i = 0; i < 4; i++) {
     1669            paramMask->data.U8[i] = 0;
     1670        }
     1671        for (int i = 4; i < paramMask->n; i++) {
     1672            paramMask->data.U8[i] = 1;
     1673        }
     1674    }
     1675
     1676    // Set the parameter range checks
     1677    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
     1678    psVector *beta_lim = NULL;
     1679    psVector *params_min = NULL;
     1680    psVector *params_max = NULL;
     1681
     1682    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
     1683    //           in the SDRS, I define a new psMinimization which will take these in
     1684    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
     1685    modelLimits (&beta_lim, &params_min, &params_max);
     1686    for (int i = 0; i < params->n; i++) {
     1687        covar->data.F64[0][i] = beta_lim->data.F32[i];
     1688        covar->data.F64[1][i] = params_min->data.F32[i];
     1689        covar->data.F64[2][i] = params_max->data.F32[i];
     1690    }
     1691
     1692    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
     1693    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     1694    for (int i = 0; i < dparams->n; i++) {
     1695        if ((paramMask != NULL) && paramMask->data.U8[i])
     1696            continue;
     1697        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     1698    }
     1699
     1700    // save the resulting chisq, nDOF, nIter
     1701    model->chisq = myMin->value;
     1702    model->nIter = myMin->iter;
     1703    model->nDOF  = y->n - nParams;
     1704
     1705    // get the Gauss-Newton distance for fixed model parameters
     1706    if (paramMask != NULL) {
     1707        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     1708        psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
    16811709        for (int i = 0; i < dparams->n; i++) {
    1682             if ((paramMask != NULL) && paramMask->data.U8[i])
     1710            if (!paramMask->data.U8[i])
    16831711                continue;
    1684             dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
    1685         }
    1686 
    1687         // save the resulting chisq, nDOF, nIter
    1688         model->chisq = myMin->value;
    1689         model->nIter = myMin->iter;
    1690         model->nDOF  = y->n - nParams;
    1691 
    1692         // get the Gauss-Newton distance for fixed model parameters
    1693         if (paramMask != NULL) {
    1694             psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
    1695             psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc);
    1696             for (int i = 0; i < dparams->n; i++) {
    1697                 if (!paramMask->data.U8[i])
    1698                     continue;
    1699                 dparams->data.F32[i] = delta->data.F64[i];
    1700             }
    1701         }
    1702 
    1703         // set the model success or failure status
    1704         if (!fitStatus) {
    1705             model->status = PM_MODEL_NONCONVERGE;
    1706         } else {
    1707             model->status = PM_MODEL_SUCCESS;
    1708         }
    1709 
    1710         // models can go insane: reject these
    1711         onPic &= (params->data.F32[2] >= source->pixels->col0);
    1712         onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    1713         onPic &= (params->data.F32[3] >= source->pixels->row0);
    1714         onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
    1715         if (!onPic) {
    1716             model->status = PM_MODEL_OFFIMAGE;
    1717         }
    1718 
    1719         source->mode |= PM_SOURCE_FITTED;
    1720 
    1721         psFree(x);
    1722         psFree(y);
    1723         psFree(yErr);
    1724         psFree(myMin);
    1725         psFree(covar);
    1726         psFree(paramMask);
    1727 
    1728         rc = (onPic && fitStatus);
    1729         psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1730         return(rc);
    1731     }
    1732 
    1733     bool p_pmSourceAddOrSubModel(psImage *image,
    1734                                  psImage *mask,
    1735                                  pmModel *model,
    1736                                  bool center,
    1737                                  bool sky,
    1738                                  bool add
    1739                                     ) {
    1740         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1741 
    1742         PS_ASSERT_PTR_NON_NULL(model, false);
    1743         PS_ASSERT_IMAGE_NON_NULL(image, false);
    1744         PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
    1745 
    1746         psVector *x = psVectorAlloc(2, PS_TYPE_F32);
    1747         psVector *params = model->params;
    1748         pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    1749         psS32 imageCol;
    1750         psS32 imageRow;
    1751         psF32 skyValue = params->data.F32[0];
    1752         psF32 pixelValue;
    1753 
    1754         for (psS32 i = 0; i < image->numRows; i++) {
    1755             for (psS32 j = 0; j < image->numCols; j++) {
    1756                 if ((mask != NULL) && mask->data.U8[i][j])
    1757                     continue;
    1758 
    1759                 // XXX: Should you be adding the pixels for the entire subImage,
    1760                 // or a radius of pixels around it?
    1761 
    1762                 // Convert i/j to imace coord space:
    1763                 // XXX: Make sure you have col/row order correct.
    1764                 // XXX EAM : 'center' option changes this
    1765                 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
    1766                 if (center) {
    1767                     imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
    1768                     imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
    1769                 } else {
    1770                     imageCol = j + image->col0;
    1771                     imageRow = i + image->row0;
    1772                 }
    1773 
    1774                 x->data.F32[0] = (float) imageCol;
    1775                 x->data.F32[1] = (float) imageRow;
    1776 
    1777                 // set the appropriate pixel value for this coordinate
    1778                 if (sky) {
    1779                     pixelValue = modelFunc (NULL, params, x);
    1780                 } else {
    1781                     pixelValue = modelFunc (NULL, params, x) - skyValue;
    1782                 }
    1783 
    1784 
    1785                 // add or subtract the value
    1786                 if (add
    1787                    ) {
    1788                     image->data.F32[i][j] += pixelValue;
    1789                 }
    1790                 else {
    1791                     image->data.F32[i][j] -= pixelValue;
    1792                 }
    1793             }
    1794         }
    1795         psFree(x);
    1796         psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
    1797         return(true);
    1798     }
    1799 
    1800 
    1801 
    1802     /******************************************************************************
    1803      *****************************************************************************/
    1804     bool pmSourceAddModel(psImage *image,
    1805                           psImage *mask,
    1806                           pmModel *model,
    1807                           bool center,
    1808                           bool sky) {
    1809         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1810         psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
    1811         psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1812         return(rc);
    1813     }
    1814 
    1815     /******************************************************************************
    1816      *****************************************************************************/
    1817     bool pmSourceSubModel(psImage *image,
    1818                           psImage *mask,
    1819                           pmModel *model,
    1820                           bool center,
    1821                           bool sky) {
    1822         psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    1823         psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
    1824         psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
    1825         return(rc);
    1826     }
    1827 
    1828     bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) {
    1829 
    1830         float obsSum = 0;
    1831         float fitSum = 0;
    1832         float sky = model->params->data.F32[0];
    1833 
    1834         pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
    1835         fitSum = modelFluxFunc (model->params);
    1836 
    1837         for (int ix = 0; ix < image->numCols; ix++) {
    1838             for (int iy = 0; iy < image->numRows; iy++) {
    1839                 if (mask->data.U8[iy][ix])
    1840                     continue;
    1841                 obsSum += image->data.F32[iy][ix] - sky;
    1842             }
    1843         }
    1844         if (obsSum <= 0)
    1845             return false;
    1846         if (fitSum <= 0)
    1847             return false;
    1848 
    1849         *fitMag = -2.5*log10(fitSum);
    1850         *obsMag = -2.5*log10(obsSum);
    1851         return (true);
    1852     }
    1853 
     1712            dparams->data.F32[i] = delta->data.F64[i];
     1713        }
     1714    }
     1715
     1716    // set the model success or failure status
     1717    if (!fitStatus) {
     1718        model->status = PM_MODEL_NONCONVERGE;
     1719    } else {
     1720        model->status = PM_MODEL_SUCCESS;
     1721    }
     1722
     1723    // models can go insane: reject these
     1724    onPic &= (params->data.F32[2] >= source->pixels->col0);
     1725    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
     1726    onPic &= (params->data.F32[3] >= source->pixels->row0);
     1727    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     1728    if (!onPic) {
     1729        model->status = PM_MODEL_OFFIMAGE;
     1730    }
     1731
     1732    source->mode |= PM_SOURCE_FITTED;
     1733
     1734    psFree(x);
     1735    psFree(y);
     1736    psFree(yErr);
     1737    psFree(myMin);
     1738    psFree(covar);
     1739    psFree(paramMask);
     1740
     1741    rc = (onPic && fitStatus);
     1742    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     1743    return(rc);
     1744}
     1745
     1746bool p_pmSourceAddOrSubModel(psImage *image,
     1747                             psImage *mask,
     1748                             pmModel *model,
     1749                             bool center,
     1750                             bool sky,
     1751                             bool add
     1752                                )
     1753{
     1754    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1755
     1756    PS_ASSERT_PTR_NON_NULL(model, false);
     1757    PS_ASSERT_IMAGE_NON_NULL(image, false);
     1758    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
     1759
     1760    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
     1761    psVector *params = model->params;
     1762    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     1763    psS32 imageCol;
     1764    psS32 imageRow;
     1765    psF32 skyValue = params->data.F32[0];
     1766    psF32 pixelValue;
     1767
     1768    for (psS32 i = 0; i < image->numRows; i++) {
     1769        for (psS32 j = 0; j < image->numCols; j++) {
     1770            if ((mask != NULL) && mask->data.U8[i][j])
     1771                continue;
     1772
     1773            // XXX: Should you be adding the pixels for the entire subImage,
     1774            // or a radius of pixels around it?
     1775
     1776            // Convert i/j to imace coord space:
     1777            // XXX: Make sure you have col/row order correct.
     1778            // XXX EAM : 'center' option changes this
     1779            // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2]
     1780            if (center) {
     1781                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
     1782                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
     1783            } else {
     1784                imageCol = j + image->col0;
     1785                imageRow = i + image->row0;
     1786            }
     1787
     1788            x->data.F32[0] = (float) imageCol;
     1789            x->data.F32[1] = (float) imageRow;
     1790
     1791            // set the appropriate pixel value for this coordinate
     1792            if (sky) {
     1793                pixelValue = modelFunc (NULL, params, x);
     1794            } else {
     1795                pixelValue = modelFunc (NULL, params, x) - skyValue;
     1796            }
     1797
     1798
     1799            // add or subtract the value
     1800            if (add
     1801               ) {
     1802                image->data.F32[i][j] += pixelValue;
     1803            }
     1804            else {
     1805                image->data.F32[i][j] -= pixelValue;
     1806            }
     1807        }
     1808    }
     1809    psFree(x);
     1810    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     1811    return(true);
     1812}
     1813
     1814
     1815
     1816/******************************************************************************
     1817 *****************************************************************************/
     1818bool pmSourceAddModel(psImage *image,
     1819                      psImage *mask,
     1820                      pmModel *model,
     1821                      bool center,
     1822                      bool sky)
     1823{
     1824    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1825    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true);
     1826    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     1827    return(rc);
     1828}
     1829
     1830/******************************************************************************
     1831 *****************************************************************************/
     1832bool pmSourceSubModel(psImage *image,
     1833                      psImage *mask,
     1834                      pmModel *model,
     1835                      bool center,
     1836                      bool sky)
     1837{
     1838    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     1839    psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false);
     1840    psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc);
     1841    return(rc);
     1842}
     1843
     1844bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask)
     1845{
     1846
     1847    float obsSum = 0;
     1848    float fitSum = 0;
     1849    float sky = model->params->data.F32[0];
     1850
     1851    pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type);
     1852    fitSum = modelFluxFunc (model->params);
     1853
     1854    for (int ix = 0; ix < image->numCols; ix++) {
     1855        for (int iy = 0; iy < image->numRows; iy++) {
     1856            if (mask->data.U8[iy][ix])
     1857                continue;
     1858            obsSum += image->data.F32[iy][ix] - sky;
     1859        }
     1860    }
     1861    if (obsSum <= 0)
     1862        return false;
     1863    if (fitSum <= 0)
     1864        return false;
     1865
     1866    *fitMag = -2.5*log10(fitSum);
     1867    *obsMag = -2.5*log10(obsSum);
     1868    return (true);
     1869}
     1870
Note: See TracChangeset for help on using the changeset viewer.