Changeset 36307

Timestamp:

Nov 23, 2013, 12:09:42 PM (13 years ago)

Author:

eugene

Message:

add gauss to 2d cache smooth; some test code in pmPCM_Minimize

Location:

branches/eam_branches/ipp-20130904

Files:

• psLib/src/imageops/psImageConvolve.h (modified) (1 diff)
• psLib/src/imageops/psImageConvolve2dCache.c (modified) (1 diff)
• psModules/src/objects/pmPCM_MinimizeChisq.c (modified) (1 diff)
• psModules/src/objects/pmPCMdata.c (modified) (1 diff)
• psphot/src/psphotExtendedSourceFits.c (modified) (1 diff)
• psphot/test/tap_psphot_galaxygrid.pro (modified) (6 diffs)

branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h

@@ -328 +328 @@
 psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma);
 bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa);
+bool psImageSmooth2dCacheKernel_Gauss (psImageSmooth2dCacheData *smdata, float sigma);
 bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata);
 

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

@@ -125 +125 @@
 }
 
+// generate a 2D smoothing kernel for supplied sigma & kappa (PS1_V1 profile).  
+bool psImageSmooth2dCacheKernel_Gauss (psImageSmooth2dCacheData *smdata, float sigma) {
+
+    // check for NULL structure elements?
+    int Ns = (int)(smdata->Nsigma * sigma);
+    Ns = PS_MAX (3, PS_MIN (Ns, 10));
+    smdata->Ns = Ns;
+
+    int Ns2 = Ns * Ns;
+
+    // we are going to use a hard-wired set of radial points
+    smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX);
+
+    float sum = 0.0;
+    for (int i = 0; i < NRAD_MAX; i++) {
+        if (radii2[i] > Ns2) {
+            smdata->radflux[i] = 0.0;
+            continue;
+        }
+        float z = 0.5 * radii2[i] / PS_SQR(sigma);
+        smdata->radflux[i] = exp(-z);
+        sum += radiiN[i] * smdata->radflux[i];
+    }
+    for (int i = 0; i < NRAD_MAX; i++) {
+        smdata->radflux[i] = smdata->radflux[i] / sum;
+    }
+
+    return true;
+}
+
 // we can use the same DATA structure on multiple images of the same size
 bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata)

branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -134 +134 @@
             psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*pcm->nPar, min->maxTol*pcm->nPar);
         }
+
+        fprintf (stderr, "%d : \n", min->iter);
+        for (int ti = 0; ti < params->n; ti++) {
+          fprintf (stderr, "%f  ", params->data.F32[ti]);
+        }
+        fprintf (stderr, "\n");
 
         char key[10]; // used for interactive responses

branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c

@@ -357 +357 @@
 
     if (modelPSF->type == modelType_GAUSS) {
-        // psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
+        psImageSmooth2dCacheKernel_Gauss (smdata, *sigma);
     }
     if (modelPSF->type == modelType_PS1_V1) {

branches/eam_branches/ipp-20130904/psphot/src/psphotExtendedSourceFits.c

@@ -315 +315 @@
 
 // XXX TEST 
-            if (!isInteractive) {
+            if (false && !isInteractive) {
                 if (!psThreadJobAddPending(job)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");

branches/eam_branches/ipp-20130904/psphot/test/tap_psphot_galaxygrid.pro

@@ -30 +30 @@
 $RefOptions = $RefOptions -nx 3000 -ny 3000
 
+if (not($?PSFMODEL))
+  $PSFMODEL = PS1_V1
+end
+
+macro reset.options
 # options for the simulated images (using the refimage for the stars)
 $FakeOptions = $BaseOptions
 $FakeOptions = $FakeOptions -exptime 30.0
 # $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
-$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1
+$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_$PSFMODEL
 $FakeOptions = $FakeOptions -nx 3000 -ny 3000
   
@@ -70 +75 @@
 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
 $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+end
+
+if (not($?FakeConfig)) reset.options
 
 list fwhm 
@@ -317 +325 @@
 end
 
+macro grid.plot.stars
+  if ($0 != 2)
+    echo "USAGE: grid.plot.stars (version)"
+    break
+  end
+  # things to examine: theta, Rmajor, AR
+
+  # go.grid.check.devexp
+
+  # check on position
+  set dX = Xot_s - int(Xin_s) - 0.5
+  set dY = Yot_s - int(Yin_s) - 0.5
+  set ARin = rin_s / Rin_s
+  set dR = Rin_s - Rot_s
+
+  create n 0 dR[]
+  set dRf = dR / Rin_s
+  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
+  label -x sequence -y "1 - R_out| / R_in|"
+  resize 700 320
+
+  # check on A.Ratio
+  set ARot = rot_s / Rot_s
+  # lim ARin ARot; clear; box; plot ARin ARot
+  set fAR = ARot / ARin
+  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
+  label -x sequence -y "AR_out| / AR_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dM = Mot_s - Min_s
+  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM    
+  label -x sequence -y "M_out| - M_in|"
+  resize 700 320
+
+  lim -n 4$1 n -1.5 1.5; clear; box; plot n dX
+  label -x sequence -y "X_out| - X_in|"
+  resize 700 320
+
+  lim -n 5$1 n -1.5 1.5; clear; box; plot n dY
+  label -x sequence -y "Y_out| - Y_in|"
+  resize 700 320
+end
+
 macro grid.mkexp.sersic
   if ($0 != 2)
@@ -550 +602 @@
   $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T
   $psphotConfig = $psphotConfig -D  PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel
-  $psphotConfig = $psphotConfig -D  PSPHOT:PSF_MODEL PS_MODEL_PS1_V1
+  $psphotConfig = $psphotConfig -D  PSPHOT:PSF_MODEL PS_MODEL_$PSFMODEL
+
+  $psphotConfig = $psphotConfig -Db PSPHOT:PSF.RESIDUALS F
+  $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F
+  $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F
 
   # ppImage / psphot on the output
@@ -649 +705 @@
   concat min min_S
   concat Min Min_S
+end
+
+macro cmf.load.stars.concat
+  if ($0 != 3)
+    echo "USAGE: cmf.load.concat (dat) (cmf)"
+    break
+  end
+
+  data $1
+  read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
+
+  subset Xin     = Xin_all     if (Type == 0)
+  subset Yin     = Yin_all     if (Type == 0)
+  subset Min     = Min_all     if (Type == 0)
+  subset RmajIn  = RmajIn_all  if (Type == 0)
+  subset RminIn  = RminIn_all  if (Type == 0)
+  subset IndexIn = IndexIn_all if (Type == 0)
+  subset Tin_rad = ThetaIn_all if (Type == 0)
+  set Tin = Tin_rad * 180 / 3.14159265
+
+  data $2
+
+  break -auto off
+  read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG PSF_MAJOR PSF_MINOR PSF_THETA
+  set PSF_THETA_ALT = PSF_THETA * (PSF_THETA >= 0.0) + (PSF_THETA + 3.14159265) * (PSF_THETA < 0.0)
+  set PSF_THETA = PSF_THETA_ALT * 180 / 3.14159265
+  
+  match2d X_PSF Y_PSF Xin Yin 1.0 -index1 index1 -index2 index2
+
+  reindex Xot_m = X_PSF using index1
+  reindex Yot_m = Y_PSF using index1
+
+  reindex Xin_m = Xin using index2
+  reindex Yin_m = Yin using index2
+
+  reindex Mot_m = PSF_INST_MAG using index1
+  reindex Tot_m = PSF_THETA using index1
+
+  reindex Min_m = Min using index2
+  reindex Tin_m = Tin using index2
+
+  reindex Rot_m = PSF_MAJOR using index1
+  reindex rot_m = PSF_MINOR using index1
+
+  reindex Rin_m = RmajIn using index2
+  reindex rin_m = RminIn using index2
+  
+  $fields = X Y M T R r
+  foreach field $fields
+    foreach set in ot
+      concat $field\$set\_m $field\$set\_s
+    end
+  end
 end
 
@@ -1144 +1253 @@
   exit 0
 end
+
+# note that t1 = original confi
+# t2 = no residuals
+# t3 = const weight PSF fit
+# t4 = const weight galaxy fits
+
+macro run.radius.loop
+  local radius Nrun 
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
+    fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4 EXP_CONV
+    $Nrun ++
+  end
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4.cmf EXP
+    set dM_m = Mot_m - Min_m
+    vstat -q dM_m
+    $mag_off = $MEDIAN
+    set dR_m = Rot_m - Rin_m
+    vstat -q dR_m
+    $rad_off = $MEDIAN
+    echo $radius $mag_off $rad_off
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end