smurf: Re-structure AST.FILT_DIFF stuff to allow it to be used in sky…

…loop

applications/smurf/libsmf/Makefile.am

@@ -166,6 +166,7 @@ smf_filter_execute.c \
 smf_filter_fromkeymap.c \
 smf_filter_getlowf.c \
 smf_filter_ident.c \
+smf_filter_mapchange.c \
 smf_filter_mce.c \
 smf_filter_notch.c \
 smf_filter_r2c.c \

applications/smurf/libsmf/smf.h.source

@@ -1071,6 +1071,9 @@ void smf_filter_getlowf( smfFilter *filt, smfHead *hdr, double *period,
 
 void smf_filter_ident( smfFilter *filt, int complex, int *status );
 
+void smf_filter_mapchange( ThrWorkForce *wf, smfDIMMData *dat,
+                           double filt_diff, int *status);
+
 void smf_filter_mce( smfFilter *filt, int noinverse, int *status );
 
 void smf_filter_notch( smfFilter *filt, const double f_low[],

applications/smurf/libsmf/smf_calcmodel_ast.c

@@ -151,6 +151,9 @@
 *        than the map itself, and uses a hard edged filter rather than a
 *        Gaussian filter. The config parameter name is changed from
 *        "gaussbg" to "filt_diff".
+*     2014-8-20 (DSB):
+*        Move ast.filt_diff stuff into smf_iteratemap, so that it can be
+*        used by skyloop.
 *     {enter_further_changes_here}
 
 *  Copyright:
@@ -220,7 +223,6 @@ void smf_calcmodel_ast( ThrWorkForce *wf,
 
   /* Local Variables */
   size_t bstride;               /* bolo stride */
-  double filt_diff;             /* Size of background filter */
   int *hitsmap;                 /* Pointer to hitsmap data */
   dim_t i;                      /* Loop counter */
   dim_t idx=0;                  /* Index within subgroup */
@@ -273,15 +275,6 @@ void smf_calcmodel_ast( ThrWorkForce *wf,
   /* Allocate job data for threads. */
   job_data = astCalloc( nw, sizeof(*job_data) );
 
-  /* Parse parameters */
-
-  /* Will we apply a smoothing constraint? */
-  astMapGet0D( kmap, "FILT_DIFF", &filt_diff );
-  if( filt_diff < 0 && *status == SAI__OK ) {
-    *status = SAI__ERROR;
-    errRep( "", FUNC_NAME ": AST.FILT_DIFF cannot be < 0.", status );
-  }
-
   /* Before applying boundary conditions, removing AST signal from residuals
      etc., flag spikes using map */
 
@@ -324,75 +317,6 @@ void smf_calcmodel_ast( ThrWorkForce *wf,
   } else {
     dat->ast_skipped = 0;
 
-    /* Constrain the change in the map. We don't if this if no previous
-       map is available. */
-    if( filt_diff > 0.0 && dat->iter > 1 ) {
-      smfData *filtermap=NULL;
-      smfFilter *filt=NULL;
-      double *change = NULL;
-
-      msgOutiff( MSG__DEBUG, "","   high-pass filtering the map-change to "
-                 "remove features larger than %g arc-sec.", status, filt_diff );
-
-      /* Form the change in the map since the previous iteration,
-         remove low frequencies from the map change, and then add the
-         remaining high frequency changes back onto the previous map to
-         get the new map. We filter the map change rather than the map
-         itself since the astronomical signal wil be far weaker (compared
-         to the spurious low frequency structures introduced by the
-         iterative process) in the map change than in the map. This means
-         we probably do not need to worry about ringing round sources etc
-         caused by the filtering process. */
-
-      /* Create a smfData to hold the changes in the map. */
-      change = astMalloc( dat->msize*sizeof( *change ) );
-      filtermap = smf_create_smfData( 0, status );
-      if( *status == SAI__OK ) {
-        filtermap->isFFT = -1;
-        filtermap->dtype = SMF__DOUBLE;
-        filtermap->pntr[0] = change;
-        filtermap->ndims = 2;
-        filtermap->lbnd[0] = dat->lbnd_out[0];
-        filtermap->lbnd[1] = dat->lbnd_out[1];
-        filtermap->dims[0] = dat->ubnd_out[0]-dat->lbnd_out[0]+1;
-        filtermap->dims[1] = dat->ubnd_out[1]-dat->lbnd_out[1]+1;
-        filtermap->hdr->wcs = astClone( dat->outfset );
-
-        /* Form the map change, setting bad values to 0. */
-        for( i=0; i<dat->msize; i++ ) {
-          if( map[i] != VAL__BADD && dat->lastmap[i] != VAL__BADD ) {
-             change[i] = map[i] - dat->lastmap[i];
-          } else {
-             change[i] = 0;
-          }
-        }
-
-        /* Create and apply a sharp-edged high-pass filter to remove
-           large-scale structures from the map change image. */
-        filt = smf_create_smfFilter( filtermap, status );
-
-        smf_filter2d_edge( filt, 1.0/filt_diff, 0, status );
-        smf_filter_execute( wf, filtermap, filt, 0, 0, status );
-
-        /* Add the remining high frequencies of the filtered map-change
-           image back onto the previous map. */
-        for( i=0; i<dat->msize; i++ ) {
-          if( map[i] != VAL__BADD && dat->lastmap[i] != VAL__BADD ) {
-             map[i] = change[i] + dat->lastmap[i];
-          } else {
-             map[i] = VAL__BADD;
-          }
-        }
-
-        /* Unset pointers to avoid freeing them */
-        filtermap->pntr[0] = NULL;
-      }
-
-      change = astFree( change );
-      smf_close_file( wf, &filtermap, status );
-      filt = smf_free_smfFilter( filt, status );
-    }
-
     /* Get a mask to apply to the map. This is determined by the "Zero_..."
        parameters in the configuration KeyMap. */
     zmask = smf_get_mask( wf, SMF__AST, keymap, dat, flags, status );

applications/smurf/libsmf/smf_filter_mapchange.c

@@ -0,0 +1,159 @@
+ /*
+*+
+*  Name:
+*     smf_filter_mapchange
+
+*  Purpose:
+*     Remove low frequency changes from the current map.
+
+*  Language:
+*     Starlink ANSI C
+
+*  Type of Module:
+*     Library routine
+
+*  Invocation:
+*     smf_filter_mapchange( ThrWorkForce *wf, smfDIMMData *dat,
+*                           double filt_diff, int *status)
+
+*  Arguments:
+*     wf = ThrWorkForce * (Given)
+*        Pointer to a pool of worker threads
+*     dat = smfDIMMData * (Given)
+*        Struct of pointers to information required by model calculation
+*     filt_diff = double (Given)
+*        The size of the filter to use, in arc-sec. No filtering is done
+*        if this is zero.
+*     status = int* (Given and Returned)
+*        Pointer to global status.
+
+*  Description:
+*     Form the change in the map since the previous iteration, remove low
+*     frequencies from the map change, and then add the remaining high
+*     frequency changes back onto the previous map to get the new map. We
+*     filter the map change rather than the map itself since the
+*     astronomical signal wil be far weaker (compared to the spurious low
+*     frequency structures introduced by the iterative process) in the map
+*     change than in the map. This means we probably do not need to worry
+*     about ringing round sources etc caused by the filtering process.
+
+*  Authors:
+*     DSB: David Berry (JAC, Hawaii)
+*     {enter_new_authors_here}
+
+*  History:
+*     20-AUG-2014 (DSB):
+*        Initial Version
+*     {enter_further_changes_here}
+
+*  Copyright:
+*     Copyright (C) 2014 Science and Technology Facilities Council.
+*     All Rights Reserved.
+
+*  Licence:
+*     This program is free software; you can redistribute it and/or
+*     modify it under the terms of the GNU General Public License as
+*     published by the Free Software Foundation; either version 3 of
+*     the License, or (at your option) any later version.
+*
+*     This program is distributed in the hope that it will be
+*     useful, but WITHOUT ANY WARRANTY; without even the implied
+*     warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+*     PURPOSE. See the GNU General Public License for more details.
+*
+*     You should have received a copy of the GNU General Public
+*     License along with this program; if not, write to the Free
+*     Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+*     MA 02110-1301, USA
+
+*  Bugs:
+*     {note_any_bugs_here}
+*-
+*/
+
+/* Starlink includes */
+#include "mers.h"
+#include "ndf.h"
+#include "sae_par.h"
+#include "star/ndg.h"
+#include "prm_par.h"
+#include "par_par.h"
+#include "star/one.h"
+
+/* SMURF includes */
+#include "libsmf/smf.h"
+
+#define FUNC_NAME "smf_filter_mapchange"
+
+void smf_filter_mapchange( ThrWorkForce *wf, smfDIMMData *dat,
+                           double filt_diff, int *status) {
+
+/* Local Variables: */
+   dim_t i;
+   double *change = NULL;
+   double *map;
+   smfData *filtermap=NULL;
+   smfFilter *filt=NULL;
+
+/* Check inherited status. Also check the filter size is non-zero. */
+   if( *status != SAI__OK || filt_diff <= 0.0 ) return;
+
+/* Store a pointer to the current map. */
+   map = dat->map;
+
+/* Tell the user what is happening. */
+   msgOutiff( MSG__DEBUG, "","   high-pass filtering the map-change to "
+              "remove features larger than %g arc-sec.", status, filt_diff );
+
+/* Create a smfData to hold the changes in the map. */
+   change = astMalloc( dat->msize*sizeof( *change ) );
+   filtermap = smf_create_smfData( 0, status );
+   if( *status == SAI__OK ) {
+      filtermap->isFFT = -1;
+      filtermap->dtype = SMF__DOUBLE;
+      filtermap->pntr[0] = change;
+      filtermap->ndims = 2;
+      filtermap->lbnd[0] = dat->lbnd_out[0];
+      filtermap->lbnd[1] = dat->lbnd_out[1];
+      filtermap->dims[0] = dat->ubnd_out[0]-dat->lbnd_out[0]+1;
+      filtermap->dims[1] = dat->ubnd_out[1]-dat->lbnd_out[1]+1;
+      filtermap->hdr->wcs = astClone( dat->outfset );
+
+/* Form the map change, setting bad values to 0. */
+      for( i = 0; i < dat->msize; i++ ) {
+         if( map[ i ] != VAL__BADD && dat->lastmap[ i ] != VAL__BADD ) {
+            change[ i ] = map[ i ] - dat->lastmap[ i ];
+         } else {
+            change[ i ] = 0;
+         }
+      }
+
+/* Create and apply a sharp-edged high-pass filter to remove large-scale
+   structures from the map change image. */
+      filt = smf_create_smfFilter( filtermap, status );
+
+      smf_filter2d_edge( filt, 1.0/filt_diff, 0, status );
+      smf_filter_execute( wf, filtermap, filt, 0, 0, status );
+
+/* Add the remining high frequencies of the filtered map-change image back
+   onto the previous map. */
+      for( i = 0; i < dat->msize; i++ ) {
+         if( map[ i ] != VAL__BADD && dat->lastmap[ i ] != VAL__BADD ) {
+            map[ i ] = change[ i ] + dat->lastmap[ i ];
+         } else {
+            map[ i ] = VAL__BADD;
+         }
+      }
+
+/* Unset pointers to avoid freeing them */
+      filtermap->pntr[0] = NULL;
+   }
+
+/* Free resources. */
+   change = astFree( change );
+   smf_close_file( wf, &filtermap, status );
+   filt = smf_free_smfFilter( filt, status );
+}
+
+
+

applications/smurf/libsmf/smf_iteratemap.c

@@ -601,6 +601,7 @@ void smf_iteratemap( ThrWorkForce *wf, const Grp *igrp, const Grp *iterrootgrp,
                      int *status ) {
 
   /* Local Variables */
+  float ast_filt_diff;          /* Size of map-change filter */
   int ast_skip;                 /* Number of iterations with no AST model */
   int bolomap=0;                /* If set, produce single bolo maps */
   size_t bstride;               /* Bolometer stride */
@@ -703,7 +704,7 @@ void smf_iteratemap( ThrWorkForce *wf, const Grp *igrp, const Grp *iterrootgrp,
   char name[GRP__SZNAM+1];      /* Buffer for storing exported model names */
   dim_t nbolo;                  /* Number of bolometers */
   size_t ncontchunks=0;         /* Number continuous chunks outside iter loop*/
-  int nhitslim;                 /* Min number of hits allowed in a map pixel */
+  int nhitslim=0;               /* Min number of hits allowed in a map pixel */
   int nm=0;                     /* Signed int version of nmodels */
   dim_t nmodels=0;              /* Number of model components / iteration */
   int noidone;                  /* Has the NOI model been calculated yet? */
@@ -722,7 +723,7 @@ void smf_iteratemap( ThrWorkForce *wf, const Grp *igrp, const Grp *iterrootgrp,
   smf_qual_t *qua_data=NULL;    /* Pointer to DATA component of qua */
   smfGroup *quagroup=NULL;      /* smfGroup of quality model files */
   int quit=0;                   /* flag indicates when to quit */
-  int rate_limited;             /* Was the MAPTOL_RATE limit hit? */
+  int rate_limited=0;           /* Was the MAPTOL_RATE limit hit? */
   int rebinflags;               /* Flags to control rebinning */
   smfArray **res=NULL;          /* Residual signal */
   double *res_data=NULL;        /* Pointer to DATA component of res */
@@ -883,8 +884,13 @@ void smf_iteratemap( ThrWorkForce *wf, const Grp *igrp, const Grp *iterrootgrp,
 
       /* A negative AST.SKIP value over-rides NUMITER. */
       ast_skip = 0;
+      ast_filt_diff = 0.0;
       if( astMapGet0A( keymap, "AST", &kmap ) ) {
          astMapGet0I( kmap, "SKIP", &ast_skip );
+
+         /* Remove low spatial frequencies in map-change? If non-zero,
+            ast.filt_diff gives the filter size, in arc-secs. */
+         astMapGet0F( kmap, "FILT_DIFF", &ast_filt_diff );
          kmap = astAnnul( kmap );
       }
 
@@ -2367,6 +2373,25 @@ void smf_iteratemap( ThrWorkForce *wf, const Grp *igrp, const Grp *iterrootgrp,
                      ": ** %f s rebinning map",
                      status, smf_timerupdate(&tv1,&tv2,status) );
 
+          /* If required, modify the map to remove low frequencies changes
+             between the new map and the old map. */
+          if( ast_filt_diff > 0.0 ) {
+
+             /* Can only do this if we have a map from the previous
+                iteration (i.e. this is not the ifrst iter). */
+             if( importsky || iter > 1 ) {
+
+                /* Do not do it on the final iteration since we want the
+                   final map to contain all the residual flux. But when
+                   run from skyloop (as indicated by numiter being 1),
+                   we always run it (skyloop will clear ast.filt_diff on
+                   the last iteration). */
+                if( quit == -1 || numiter == 1 ) {
+                   smf_filter_mapchange( wf, &dat, ast_filt_diff, status );
+                }
+             }
+          }
+
           /* If required, subtract an external error map from the map
              created above. */
           if( emapdata && !dat.ast_skipped ) {

applications/smurf/scripts/skyloop.py

@@ -502,6 +502,13 @@ def cleanup():
       cleanup()
       sys.exit()
 
+#  See if low frequency changes are to be removed from the map on each
+#  iteration.
+   ast_filt_diff = float( invoke( "$KAPPA_DIR/configecho "
+                               "name=ast.filt_diff config={0} "
+                               "defaults=$SMURF_DIR/smurf_makemap.def "
+                               "select=\"\'450=0,850=1\'\" defval=0.0".format(config)))
+
 #  The first invocation of makemap will create NDFs holding cleaned
 #  time-series data, EXT, LUT and NOI model values. The NDFs are created
 #  with hard-wired names and put in the current working directory. For
@@ -766,12 +773,18 @@ def cleanup():
 
 #  Also, if this is the last iteration, create a modified configuration file
 #  that supresses masking (unless the xxx.zero_notlast value in the
-#  supplied config indicates otjerwise).
+#  supplied config indicates otherwise).
             for model in ["ast", "com", "flt"]:
                if zero_notlast[model] != 0:
                   add["ast.zero_notlast"] = 1
                   newcon = 1
 
+#  Also, if this is the last iteration, do not remove low frequency
+#  changes from the map.
+            if ast_filt_diff != 0.0:
+               add["ast.filt_diff"] = 0.0
+               newcon = 1
+
 #  Also override the normal values for parameters that have a
 #  corresponding "_last" value.
             if filt_edge_largescale_last != -1.0: