Changeset 42584

Timestamp:

Jan 13, 2024, 2:39:37 PM (3 years ago)

Author:

eugene

Message:

tweaks to cdensify, still not ideal

File:

• branches/eam_branches/ipp-20230313/Ohana/src/opihi/cmd.astro/cdensify.c (modified) (8 diffs)

branches/eam_branches/ipp-20230313/Ohana/src/opihi/cmd.astro/cdensify.c

@@ -1 +1 @@
 # include "data.h"
+# define dCOS(A)   ((double) cos ((double)RAD_DEG*A))
 
 # define CHECKVAL(ARG) if (!isfinite(ARG)) { gprint (GP_ERR, "illegal value for %s: %f\n", #ARG, ARG); return (FALSE); }
@@ -116 +117 @@
   PutCoords (&newcoords, &bf[0].header);
   
+  // use the mask to prevent double-counting
+  ALLOCATE_PTR (mask, char, Nx*Ny); 
+
   float scalescale = scale*scale;
   float scale2 = (scale + 1.0) * (scale + 1.0);
@@ -142 +146 @@
     coords.crval1 = rn;
     coords.crval2 = *d;
+    
+    // XXX do not oversample by more than a factor of 10
+    float dXn = dX * MAX(fabs(dCOS(*d)), 0.1);
 
     float F = 1.0;
     if (vv) { F = isFloatScale ? Fs[i] : Is[i]; }
+
+    // reset the mask so we do not double-count
+    memset (mask, 0, Nx*Ny);
 
     switch (PSFTYPE) {
@@ -162 +172 @@
         break;
       case IS_SQUARE:
-        for (ix = -scale; ix <= scale; ix += dX) {
+        for (ix = -scale; ix <= scale; ix += dXn) {
           for (iy = -scale; iy <= scale; iy += dY) {
             double rp, dp;
@@ -184 +194 @@
         break;
       case IS_CIRCLE:
-        for (ix = -scale; ix <= scale; ix += dX) {
-          for (iy = -scale; iy <= scale; iy += dY) {
+        for (ix = -scale - dXn; ix <= scale + dXn; ix += dXn) {
+          for (iy = -scale - dY; iy <= scale + dY; iy += dY) {
             float r2 = ix*ix + iy*iy;
             double rp, dp;
@@ -199 +209 @@
             if (Xb < 0) continue;
             if (Yb < 0) continue;
+            if (mask[Xb + Yb*Nx]) continue;
+            mask[Xb + Yb*Nx] = 1;
             if (vv) {
               val[Xb + Yb*Nx] += Normalize ? fCircle*F : F;
@@ -208 +220 @@
         break;
       case IS_GAUSS:
-        for (ix = -3.0*scale; ix <= 3.0*scale; ix += dX) {
+        for (ix = -3.0*scale; ix <= 3.0*scale; ix += dXn) {
           for (iy = -3.0*scale; iy <= 3.0*scale; iy += dY) {
             float r2 = ix*ix + iy*iy;
@@ -232 +244 @@
     }
   }
+
+  free (mask);
   return (TRUE);
 }