Changeset 5048 for trunk/psphot/src/pmObjects_EAM.c

Timestamp:

Sep 13, 2005, 3:35:20 PM (21 years ago)

Author:

eugene

Message:

substantial work to match with psLib v.7

File:

• trunk/psphot/src/pmObjects_EAM.c (modified) (21 diffs)

trunk/psphot/src/pmObjects_EAM.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-08 04:32:40 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-14 01:35:20 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -603 +603 @@
      members.
 *****************************************************************************/
-pmSource *pmSourceLocalSky(const psImage *image,
-                           const pmPeak *peak,
-                           psStatsOptions statsOptions,
-                           psF32 innerRadius,
-                           psF32 outerRadius)
-{
-    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL);
-    PS_ASSERT_PTR_NON_NULL(peak, NULL);
-    PS_FLOAT_COMPARE(0.0, innerRadius, NULL);
-    PS_FLOAT_COMPARE(innerRadius, outerRadius, NULL);
-    psS32 innerRadiusS32 = (psS32) innerRadius;
-    psS32 outerRadiusS32 = (psS32) outerRadius;
-
-    //
-    // We define variables for code readability.
-    //
-    // XXX: Since the peak->xy coords are in image, not subImage coords,
-    // these variables should be renamed for clarity (imageCenterRow, etc).
-    //
-    // peak->x,y is guaranteed to be on image
-    psS32 SubImageCenterRow = peak->y;
-    psS32 SubImageCenterCol = peak->x;
-
-    // XXX EAM : I added this code to stay on the image. So did George
-    // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1)
-    // XXX EAM : dropped the Annulus width (not needed)
-    // XXX EAM : dropped off-boundary tests (not needed)
-    psS32 SubImageStartRow  = PS_MAX (0, SubImageCenterRow - outerRadiusS32);
-    psS32 SubImageEndRow    = PS_MIN (image->numRows, SubImageCenterRow + outerRadiusS32 + 1);
-    psS32 SubImageStartCol  = PS_MAX (0, SubImageCenterCol - outerRadiusS32);
-    psS32 SubImageEndCol    = PS_MIN (image->numCols, SubImageCenterCol + outerRadiusS32 + 1);
-    //
-    // Grab a subimage of the original image of size (2 * outerRadius).
-    //
-    // XXX EAM : merged psImageSubset & psRegionSet
-    // XXX EAM : cast for image is needed because of const (above)
-    psImage *subImage = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol,
-                                                         SubImageEndCol,
-                                                         SubImageStartRow,
-                                                         SubImageEndRow));
-
-    psImage *subImageMask = psImageAlloc(subImage->numCols,
-                                         subImage->numRows,
-                                         PS_TYPE_U8);
-    // XXX EAM : for consistency, mask needs col0,row0 set to match image
-    // XXX EAM : CONFLICT between psLib and code requirement. FIX PSLIB!!!
-    // subImageMask->col0 = subImage->col0;
-    // subImageMask->row0 = subImage->row0;
-    psAbort ("pmObjects", "must fix this error before psphot will work!!!\n");
-
-    //
-    // Loop through the subimage mask & initialize PSPHOT_MASK_INVALID
-    //
-    for (psS32 row = 0 ; row < subImageMask->numRows; row++) {
-        for (psS32 col = 0 ; col < subImageMask->numCols; col++) {
-            subImageMask->data.U8[row][col] = PSPHOT_MASK_INVALID;
-        }
-    }
-
-    //
-    // Loop through the subimage, mask off pixels in the inner square.
-    // XXX this uses a static mask value of 0
-    // XXX EAM : this was wrong: it masked the wrong pixels at the edge of the image
-    // XXX EAM : this masks the pixels in the center region
-    // XXX EAM : should we be using psImageMaskRegion ???
-    psS32 StartRow  = PS_MAX (0, SubImageCenterRow - subImageMask->row0 - innerRadiusS32);
-    psS32 EndRow    = PS_MIN (subImageMask->numRows, SubImageCenterRow - subImageMask->row0 + innerRadiusS32 + 1);
-    psS32 StartCol  = PS_MAX (0, SubImageCenterCol - subImageMask->col0 - innerRadiusS32);
-    psS32 EndCol    = PS_MIN (subImageMask->numCols, SubImageCenterCol - subImageMask->col0 + innerRadiusS32 + 1);
-    for (psS32 row = StartRow; row < EndRow; row++) {
-        for (psS32 col = StartCol; col < EndCol; col++) {
-            subImageMask->data.U8[row][col] = PSPHOT_MASK_CLEAR;
-        }
-    }
-
-    // XXX EAM : make this trace information??
-    //    for (psS32 row = 0 ; row < subImage->numRows; row++) {
-    //        for (psS32 col = 0 ; col < subImage->numCols; col++) {
-    //            printf("(%d) ", subImageMask->data.U8[row][col]);
-    //        }
-    //        printf("\n");
-    //    }
-
-    //
-    // Allocate the myStats structure, then call psImageStats(), which will
-    // calculate the specified statistic.
-    //
+bool pmSourceLocalSky (pmSource *source,
+                           psStatsOptions statsOptions,
+                           psF32 Radius)
+{
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    // XXX EAM : psRegionXXX funcs should take value not ptr
+    srcRegion = psRegionForSquare (peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage (mask, &srcRegion);
+
+    psImageMaskRegion (mask, &srcRegion, "OR", PSPHOT_MASK_MARKED);
     psStats *myStats = psStatsAlloc(statsOptions);
-    myStats = psImageStats(myStats, subImage, subImageMask, 1);
-
-    //
-    // Create the output mySource, and set appropriate members.
-    //
-    pmSource *mySource = pmSourceAlloc();
-    mySource->peak = (pmPeak *) peak;
-    mySource->moments = pmMomentsAlloc();
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion (mask, &srcRegion, "AND", ~PSPHOT_MASK_MARKED);
+
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
-    mySource->moments->Sky = (psF32) tmpF64;
-    mySource->pixels = subImage;
-    mySource->mask = subImageMask;
-
-    //
-    // Free things.  XXX: This should be static memory.
-    //
-    psFree(myStats);
-
-    return(mySource);
-}
-
-/******************************************************************************
-bool checkRadius(*peak, radius, x, y): private function which simply
-determines if the (x, y) point is within the radius of the specified peak.
- 
-XXX: macro this for performance.
-*****************************************************************************/
-static bool checkRadius(pmPeak *peak,
-                        psF32 radius,
-                        psS32 x,
-                        psS32 y)
-{
-    if (PS_SQR(radius) >= (psF32) (PS_SQR(x - peak->x) + PS_SQR(y - peak->y))) {
-        return(true);
-    }
-
-    return(false);
+    psFree (myStats);
+
+    if (isnan(tmpF64)) return (false);
+    source->moments = pmMomentsAlloc();
+    source->moments->Sky = (psF32) tmpF64;
+    return (true);
 }
 
@@ -762 +662 @@
     pmSource->peak
     pmSource->pixels
+    pmSource->noise
+    pmSource->mask
  
 XXX: The peak calculations are done in image coords, not subImage coords.
  
-XXX: mask values?
-*****************************************************************************/
+XXX EAM : this version clips input pixels on S/N
+XXX EAM : this version returns false for several reasons
+*****************************************************************************/
+# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
+
 bool pmSourceMoments(pmSource *source,
                      psF32 radius)
 {
-    PS_ASSERT_PTR_NON_NULL(source, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->peak, NULL);
-    PS_ASSERT_PTR_NON_NULL(source->pixels, NULL);
+    // PS_PTR_CHECK_NULL(source, NULL);
+    // PS_PTR_CHECK_NULL(source->peak, NULL);
+    // PS_PTR_CHECK_NULL(source->pixels, NULL);
+    // PS_PTR_CHECK_NULL(source->noise, NULL);
     PS_FLOAT_COMPARE(0.0, radius, NULL);
 
@@ -791 +697 @@
     // XY  = SUM (x - xc)*(y - yc)*(z - sky)
     //
-    psF32 Sum = 0.0;
     psF32 peakPixel = -PS_MAX_F32;
     psS32 numPixels = 0;
+    psF32 Sum = 0.0;
     psF32 X1 = 0.0;
     psF32 Y1 = 0.0;
@@ -799 +705 @@
     psF32 Y2 = 0.0;
     psF32 XY = 0.0;
-    psF32 x = 0;
-    psF32 y = 0;
-    //
+    psF32 x  = 0;
+    psF32 y  = 0;
+    psF32 R2 = PS_SQR(radius);
+
+    psF32 xPeak = source->peak->x;
+    psF32 yPeak = source->peak->y;
+
     // XXX why do I get different results for these two methods of finding Sx?
     // XXX Sx, Sy would be better measured if we clip pixels close to sky
@@ -809 +719 @@
     // proceed with the moments calculation.  need to do two loops for a
     // numerically stable result.  first loop: get the sums.
-    //
+    // XXX EAM : mask == 0 is valid
+
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
         for (psS32 col = 0; col < source->pixels->numCols ; col++) {
-            // XXX EAM : mask == 0 is valid
             if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
-            psS32 imgColCoord = col + source->pixels->col0;
-            psS32 imgRowCoord = row + source->pixels->row0;
-            if (checkRadius(source->peak,
-                            radius,
-                            imgColCoord,
-                            imgRowCoord)) {
-                psF32 xDiff = (psF32) (imgColCoord - source->peak->x);
-                psF32 yDiff = (psF32) (imgRowCoord - source->peak->y);
-                psF32 pDiff = source->pixels->data.F32[row][col] - sky;
-
-                Sum+= pDiff;
-                X1+= xDiff * pDiff;
-                Y1+= yDiff * pDiff;
-                XY+= xDiff * yDiff * pDiff;
-
-                X2+= PS_SQR(xDiff) * pDiff;
-                Y2+= PS_SQR(yDiff) * pDiff;
-
-                if (source->pixels->data.F32[row][col] > peakPixel) {
-                    peakPixel = source->pixels->data.F32[row][col];
-                }
-                numPixels++;
-            }
+
+            psF32 xDiff = col + source->pixels->col0 - xPeak;
+            psF32 yDiff = row + source->pixels->row0 - yPeak;
+
+            // XXX EAM : calculate xDiff, yDiff up front;
+            //           radius is just a function of (xDiff, yDiff)
+            if (!VALID_RADIUS(xDiff, yDiff, R2)) continue;
+
+            psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+
+            // XXX EAM : check for valid S/N in pixel
+            // XXX EAM : should this limit be user-defined?
+            if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue;
+            
+            Sum += pDiff;
+            X1  += xDiff * pDiff;
+            Y1  += yDiff * pDiff;
+            XY  += xDiff * yDiff * pDiff;
+            
+            X2  += PS_SQR(xDiff) * pDiff;
+            Y2  += PS_SQR(yDiff) * pDiff;
+            
+            peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
+            numPixels++;
         }
     }
     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
     if ((numPixels < 3) || (Sum <= 0)) {
-        psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
-        return (false);
+      psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+      return (false);
     }
 
@@ -857 +768 @@
     y = Y1/Sum;
     if ((fabs(x) > radius) || (fabs(y) > radius)) {
-        psTrace (".psModules.pmSourceMoments", 5, 
-                 "large centroid swing; invalid peak %d, %d\n", 
-                 source->peak->x, source->peak->y);
-        return (false);
-    }
-
-    source->moments->x = x + ((psF32) source->peak->x);
-    source->moments->y = y + ((psF32) source->peak->y);
-
+      psTrace (".psModules.pmSourceMoments", 5, 
+               "large centroid swing; invalid peak %d, %d\n", 
+               source->peak->x, source->peak->y);
+      return (false);
+    }
+
+    source->moments->x = x + xPeak;
+    source->moments->y = y + yPeak;
+
+    // XXX EAM : Sxy needs to have x*y subtracted
     source->moments->Sxy = XY/Sum - x*y;
     source->moments->Sum = Sum;
@@ -919 +831 @@
 *****************************************************************************/
 
+// XXX EAM include a S/N cutoff in selecting the sources?
 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata)
 {
@@ -1123 +1036 @@
         bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY));
         if ((Nsatpix > 1) && big) {
-            tmpSrc->type |= PS_SOURCE_SATSTAR;
+            tmpSrc->type |= PM_SOURCE_SATSTAR;
             Nsatstar ++;
             continue;
@@ -1130 +1043 @@
         // saturated object (not a star, eg bleed trails, hot pixels)
         if (Nsatpix > 1) {
-            tmpSrc->type |= PS_SOURCE_SATURATED;
+            tmpSrc->type |= PM_SOURCE_SATURATED;
             Nsat ++;
             continue;
@@ -1139 +1052 @@
         // only set candidate defects if 
         if ((sigX < 0.05) || (sigY < 0.05)) {
-            tmpSrc->type |= PS_SOURCE_DEFECT;
+            tmpSrc->type |= PM_SOURCE_DEFECT;
             Ncr ++;
             continue;
@@ -1146 +1059 @@
         // likely unsaturated galaxy (too large to be stellar)
         if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) {
-            tmpSrc->type |= PS_SOURCE_GALAXY;
+            tmpSrc->type |= PM_SOURCE_GALAXY;
             Ngal ++;
             continue;
@@ -1159 +1072 @@
         psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY);
         if ((SN > PSF_SN_LIM) && (radius < 1.5)) {
-            tmpSrc->type |= PS_SOURCE_PSFSTAR;
+            tmpSrc->type |= PM_SOURCE_PSFSTAR;
             Npsf ++;
             continue;
@@ -1165 +1078 @@
 
         // random type of star
-        tmpSrc->type |= PS_SOURCE_OTHER;
+        tmpSrc->type |= PM_SOURCE_OTHER;
     }
 
@@ -1197 +1110 @@
 XXX: The circle will have a diameter of (1+radius).  This is different from
      the pmSourceSetLocal() function.
+
+XXX: this should probably use the psImageMaskCircle Function for consistency
 *****************************************************************************/
 bool pmSourceSetPixelsCircle(pmSource *source,
@@ -1537 +1452 @@
 XXX EAM : fit the specified model (not necessarily the one in source)
 *****************************************************************************/
-bool pmSourceFitModel(pmSource *source,
-                      pmModel *model,
-                      const bool PSF)
+bool pmSourceFitModel_v5(pmSource *source,
+                         pmModel *model,
+                         const bool PSF)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -1554 +1469 @@
     //           tests below could be conditions (!NULL)
 
-    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
-
     psVector *params = model->params;
     psVector *dparams = model->dparams;
     psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
 
     int nParams = PSF ? params->n - 4 : params->n;
@@ -1615 +1530 @@
     }       
 
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+    // XXX EAM : covar must be F64?
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
@@ -1621 +1537 @@
     for (int i = 0; i < dparams->n; i++) {
         if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
-        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     }
  
@@ -1640 +1556 @@
     // XXX EAM get the Gauss-Newton distance for fixed model parameters
     if (paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
         psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
         for (int i = 0; i < dparams->n; i++) {
             if (!paramMask->data.U8[i]) continue;
-            dparams->data.F32[i] = delta->data.F32[i];
+            dparams->data.F32[i] = delta->data.F64[i];
         }
         psFree (delta);
@@ -1655 +1571 @@
     psFree(covar);
     psFree(paramMask);
+
+    rc = (onPic && fitStatus);
+    return(rc);
+}
+
+// XXX EAM : new version with parameter range limits and noise enhancement
+bool pmSourceFitModel (pmSource *source,
+                       pmModel *model,
+                       const bool PSF)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->moments, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->noise, false);
+
+    // XXX EAM : is it necessary for the mask & noise to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+    psF32  Ro, ymodel;
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+
+    // XXX EAM : I need to use the sky value to constrain the noise model
+    int nParams = PSF ? params->n - 4 : params->n;
+    psF32 So = params->data.F32[0];
+
+    // find the number of valid pixels
+    // XXX EAM : this loop and the loop below could just be one pass
+    //           using the psArrayAdd and psVectorExtend functions
+    psS32 count = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                count++;
+            }
+        }
+    }
+    if (count <  nParams + 1) {
+        psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+        return(false);
+    }
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(count);
+    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
+    psS32 tmpCnt = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+                // XXX: Convert i/j to image space:
+                // XXX EAM: coord order is (x,y) == (col,row)
+                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+                x->data[tmpCnt] = (psPtr *) coord;
+                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
+
+                // compare observed flux to model flux to adjust weight
+                ymodel = modelFunc (NULL, model->params, coord);
+                
+                // this test enhances the noise based on deviation from the model flux
+                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+
+                // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 
+                //           the minimization function calculates sq(). 
+                //           should psMinimizeLMChi2 take dY^2?
+                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
+                tmpCnt++;
+            }
+        }
+    }
+
+    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
+                                                PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+        paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+        for (int i = 0; i < 4; i++) {
+            paramMask->data.U8[i] = 0;
+        }
+        for (int i = 4; i < paramMask->n; i++) {
+            paramMask->data.U8[i] = 1;
+        }
+    }  
+
+    // XXX EAM : I've added three types of parameter range checks
+    // XXX EAM : this requires my new psMinimization functions
+    pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
+    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
+    //           in the SDRS, I define a new psMinimization which will take these in
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+        covar->data.F64[0][i] = beta_lim->data.F32[i];
+        covar->data.F64[1][i] = params_min->data.F32[i];
+        covar->data.F64[2][i] = params_max->data.F32[i];
+    }
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
+        dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+ 
+    // XXX EAM: we need to do something (give an error?) if rc is false
+    // XXX EAM: psMinimizeLMChi2 does not check convergence
+
+    // XXX models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+
+    // XXX EAM: save the resulting chisq, nDOF, nIter
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+        psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
+        for (int i = 0; i < dparams->n; i++) {
+            if (!paramMask->data.U8[i]) continue;
+            dparams->data.F32[i] = delta->data.F64[i];
+        }
+    }
+
+    psFree(paramMask);
+    psFree(x);
+    psFree(y);
+    psFree(myMin);
 
     rc = (onPic && fitStatus);