From a4afb88b08d4efb3d45ae3f997e347d300eb4886 Mon Sep 17 00:00:00 2001 From: David Berry Date: Wed, 20 Aug 2014 13:41:50 +0100 Subject: [PATCH] smurf: Re-structure AST.FILT_DIFF stuff to allow it to be used in skyloop --- applications/smurf/libsmf/Makefile.am | 1 + applications/smurf/libsmf/smf.h.source | 3 + applications/smurf/libsmf/smf_calcmodel_ast.c | 82 +-------- .../smurf/libsmf/smf_filter_mapchange.c | 159 ++++++++++++++++++ applications/smurf/libsmf/smf_iteratemap.c | 29 +++- applications/smurf/scripts/skyloop.py | 15 +- 6 files changed, 207 insertions(+), 82 deletions(-) create mode 100644 applications/smurf/libsmf/smf_filter_mapchange.c diff --git a/applications/smurf/libsmf/Makefile.am b/applications/smurf/libsmf/Makefile.am index f48faf70a79..285f756081e 100644 --- a/applications/smurf/libsmf/Makefile.am +++ b/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 \ diff --git a/applications/smurf/libsmf/smf.h.source b/applications/smurf/libsmf/smf.h.source index 1d9479230ad..5845cdef4c0 100644 --- a/applications/smurf/libsmf/smf.h.source +++ b/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[], diff --git a/applications/smurf/libsmf/smf_calcmodel_ast.c b/applications/smurf/libsmf/smf_calcmodel_ast.c index 7b851f2a8b0..6b5b116270a 100644 --- a/applications/smurf/libsmf/smf_calcmodel_ast.c +++ b/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; imsize; 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; imsize; 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 ); diff --git a/applications/smurf/libsmf/smf_filter_mapchange.c b/applications/smurf/libsmf/smf_filter_mapchange.c new file mode 100644 index 00000000000..3390c10bb1d --- /dev/null +++ b/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 ); +} + + + diff --git a/applications/smurf/libsmf/smf_iteratemap.c b/applications/smurf/libsmf/smf_iteratemap.c index 3ce33cfd62d..3a023413d00 100644 --- a/applications/smurf/libsmf/smf_iteratemap.c +++ b/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 ) { diff --git a/applications/smurf/scripts/skyloop.py b/applications/smurf/scripts/skyloop.py index 829fa1e474e..e143716359a 100755 --- a/applications/smurf/scripts/skyloop.py +++ b/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: