smurf: More fixes to support mask freezing in skyloop

applications/smurf/libsmf/smf_initial_sky.c

@@ -71,6 +71,11 @@
 *     10-MAR-2014 (DSB):
 *        Update the Quality component of the supplied NDF to match the
 *        mask created by smf_calcmodel_ast.
+*     8-SEP-2014 (DSB):
+*        Changed to support freezing of masks within SKYLOOP. Initial
+*        masks are now made from the Quality array in the supplied NDF.
+*        These will be immediately over-written (within smf_calcmodel_ast)
+*        with new masks unless the mask is frozen.
 *     {enter_further_changes_here}
 
 *  Copyright:
@@ -122,6 +127,8 @@ int smf_initial_sky( ThrWorkForce *wf, AstKeyMap *keymap, smfDIMMData *dat,
    int result;                /* Returned flag */
    int there;                 /* Is there a smurf extension in the NDF? */
    int update;                /* Was NDF opened for UPDATE access? */
+   size_t i;                  /* Loop count */
+   size_t junk;               /* Unused value */
 
 /* Initialise the returned value to indicate no sky has been subtractred. */
    result = 0;
@@ -205,6 +212,37 @@ int smf_initial_sky( ThrWorkForce *wf, AstKeyMap *keymap, smfDIMMData *dat,
       if( there ) ndfXgt0i( indf1, SMURF__EXTNAME, "NUMITER", iters,
                             status );
 
+/* If the NDF has a Quality component, import it and create initial AST,
+   FLT and COM masks from it. These will often be over-ridden by new masks
+   calculated with smf_calcmodel_ast below, but will not be over-written
+   if the masks have been frozen by xxx.zero_freeze. */
+      ndfState( indf2, "Quality", &there, status );
+      if( there && dat->mapqual ) {
+         smf_qual_t *qarray = smf_qual_map( wf, indf2, "Read", NULL, &junk,
+                                            status );
+         if( *status == SAI__OK ) {
+            smf_qual_t *pq = qarray;
+            for( i = 0; i < dat->msize; i++,pq++ ) {
+               if( *pq & SMF__MAPQ_AST ) {
+                  if( !dat->ast_mask ) dat->ast_mask = astCalloc( dat->msize,
+                                                  sizeof( *(dat->ast_mask) ) );
+                  (dat->ast_mask)[ i ] = 1;
+               }
+               if( *pq & SMF__MAPQ_FLT ) {
+                  if( !dat->flt_mask ) dat->flt_mask = astCalloc( dat->msize,
+                                                  sizeof( *(dat->flt_mask) ) );
+                  (dat->flt_mask)[ i ] = 1;
+               }
+               if( *pq & SMF__MAPQ_COM ) {
+                  if( !dat->com_mask ) dat->com_mask = astCalloc( dat->msize,
+                                                  sizeof( *(dat->com_mask) ) );
+                  (dat->com_mask)[ i ] = 1;
+               }
+            }
+         }
+         qarray = astFree( qarray );
+      }
+
 /* Indicate the map arrays within the supplied smfDIMMData structure now
    contain usable values. We need to do this before calling
    smf_calcmodel_ast below so that the right mask gets used in

applications/smurf/scripts/skyloop.py

@@ -772,7 +772,7 @@ def cleanup():
                newcon = 1
 
 #  When "com_freeze_flags" invocations have been performed, freeze the
-#  COM flags (so long as com_freeze_flags > 0).  Do this for AST, COM and FLT models.
+#  COM flags (so long as com_freeze_flags > 0).
          if com_freeze_flags > 0 and iter > com_freeze_flags + 1:
             com_freeze_flags = 0
             add[ "com.freeze_flags" ] = -1
@@ -867,11 +867,26 @@ def cleanup():
                cmd += " "+extra
             invoke(cmd)
 
+#  The quality array in the new map will not be of much use since it will
+#  have been created on the basis of maps made from individual chunks, rather
+#  than the total coadded map. This will causes the following estimation
+#  of the normalised change to be wrong. So we copy the quality mask from
+#  the previous map to the new map, and use that instead (this mask was
+#  created when the previous map was read into makemap). This also helps
+#  if the mask is frozen by one of the xxx.zero_freeze config parameters.
+            if prevmap != None:
+               try:
+                  invoke("$KAPPA_DIR/setqual ndf={0} like={1}".format(newmap,prevmap) )
+               except starutil.StarUtilError as err:
+                  pass
+
 #  If required, get the mean normalised map change, and see if it has
 #  dropped below maptol. If so, we must do one further iteration to
 #  ensure that the masking is not visible in the final map.
+         invoke("$KAPPA_DIR/setbb ndf={0} bb=1".format(newmap) )
          invoke("$KAPPA_DIR/maths exp=\"'abs(ia-ib)/sqrt(va)'\" ia={0} "
                 "ib={1} out={2}".format(newmap,prevmap,mapchange))
+         invoke("$KAPPA_DIR/setbb ndf={0} bb=0".format(newmap) )
          invoke("$KAPPA_DIR/stats ndf={0} clip=\\[3,3,3\\] quiet".format(mapchange))
          meanchange = starutil.get_task_par( "mean", "stats" )
          if maptol > 0.0 and converged == False:
@@ -903,17 +918,6 @@ def cleanup():
       inputs = NDG( maps )
       invoke("$KAPPA_DIR/paste in={0} out={1} shift=\[0,0,1\]".format(inputs,itermap) )
 
-#  Since no masking is done on the final iteration, the output map will
-#  not currently contain a quality array. If the previous iteration has a
-#  quality component copy it to the last iteration, and set bad bits to
-#  zero. This is to provide a record of the final used mask.
-   if prevmap != None:
-      try:
-         invoke("$KAPPA_DIR/setqual ndf={0} like={1}".format(newmap,prevmap) )
-      except starutil.StarUtilError as err:
-         pass
-   if niter > 1 : invoke("$KAPPA_DIR/setbb ndf={0} bb=0".format(newmap) )
-
 #  Remove temporary files.
    cleanup()