debian/changelog

-snaphu (2.0.0-2) UNRELEASED; urgency=medium
+snaphu (2.0.1-1) unstable; urgency=medium
 
   [ Bas Couwenberg ]
   * Don't delete bin directory in clean target, included in upstream source.
 
   [ Antonio Valentino ]
+  * New upstream release.
   * debian/tests/control:
     - mark test as superficial using Restrictions
+  * debian/patches:
+    - refresh all patches
+  * Bump debhelper from old 11 to 12.
+  * Remove obsolete fields Name from debian/upstream/metadata.
 
- -- Bas Couwenberg <sebastic@debian.org>  Wed, 10 Jul 2019 19:47:39 +0200
+ -- Antonio Valentino <antonio.valentino@tiscali.it>  Sun, 08 Sep 2019 11:08:52 +0000
 
 snaphu (2.0.0-1) unstable; urgency=medium
 

debian/compat

-11

debian/control

@@ -4,7 +4,7 @@ Uploaders: Antonio Valentino <antonio.valentino@tiscali.it>
 Section: non-free/science
 XS-Autobuild: no
 Priority: optional
-Build-Depends: debhelper (>= 11)
+Build-Depends: debhelper-compat (= 12)
 Standards-Version: 4.4.0
 Vcs-Browser: https://salsa.debian.org/debian-gis-team/snaphu
 Vcs-Git: https://salsa.debian.org/debian-gis-team/snaphu.git

debian/patches/0002-Spelling.patch

@@ -21,10 +21,10 @@ index 161cc25..a4ed3c9 100644
  preceding one, although round-off errors in flow-to-phase conversions
  may cause minor differences
 diff --git a/src/snaphu.h b/src/snaphu.h
-index 3a9ab8d..6ec5974 100644
+index 5d43cb5..82e6d2b 100644
 --- a/src/snaphu.h
 +++ b/src/snaphu.h
-@@ -406,7 +406,7 @@
+@@ -407,7 +407,7 @@
   "The parts of this software derived from the CS2 minimum cost flow\n"\
   "solver written by A. V. Goldberg and B. Cherkassky are governed by the\n"\
   "terms of the copyright holder of that software.  Permission has been\n"\

debian/upstream/metadata

----
-Name: SNAPHU
 Other-References: https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/
 

src/snaphu.c

@@ -254,6 +254,9 @@ int Unwrap(infileT *infiles, outfileT *outfiles, paramT *params,
               /* see if next tile needs to be unwrapped */
               if(dotilemask[nexttilerow][nexttilecol]){
 
+                /* wait to make sure file i/o, threads, and OS are synched */
+                sleep(sleepinterval);
+                
                 /* fork to create new process */
                 fflush(NULL);
                 pid=fork();
@@ -317,10 +320,9 @@ int Unwrap(infileT *infiles, outfileT *outfiles, paramT *params,
                 nexttilerow++;
               }
 
-              /* wait a little while for file i/o before beginning next tile */
+              /* increment counter of running child processes */
               if(pid!=iterparams->parentpid){
                 nchildren++;
-                sleep(sleepinterval);
               }
 
             }else{
@@ -414,13 +416,14 @@ int UnwrapTile(infileT *infiles, outfileT *outfiles, paramT *params,
   long nflow, ncycle, mostflow, nflowdone;
   long candidatelistsize, candidatebagsize;
   long isource, nsource;
+  long nincreasedcostiter;
   long *nconnectedarr;
   int *nnodesperrow, *narcsperrow;
   short **flows, **mstcosts;
   float **wrappedphase, **unwrappedphase, **mag, **unwrappedest;
   incrcostT **incrcosts;
   void **costs;
-  totalcostT totalcost, oldtotalcost;
+  totalcostT totalcost, oldtotalcost, mintotalcost;
   nodeT **sourcelist;
   nodeT *source, ***apexes;
   nodeT **nodes, ground[1];
@@ -541,6 +544,9 @@ int UnwrapTile(infileT *infiles, outfileT *outfiles, paramT *params,
 	      &candidatelist,&iscandidate,&apexes,&bkts,&iincrcostfile,
 	      &incrcosts,&nodes,ground,&nnoderow,&nnodesperrow,&narcrow,
 	      &narcsperrow,nrow,ncol,&notfirstloop,&totalcost,params);
+  oldtotalcost=totalcost;
+  mintotalcost=totalcost;
+  nincreasedcostiter=0;
 
   /* regrow regions with -G parameter */
   if(params->regrowconncomps){
@@ -633,10 +639,17 @@ int UnwrapTile(infileT *infiles, outfileT *outfiles, paramT *params,
       if(notfirstloop){
 	oldtotalcost=totalcost;
 	totalcost=EvaluateTotalCost(costs,flows,nrow,ncol,NULL,params);
+        if(totalcost<mintotalcost){
+          mintotalcost=totalcost;
+        }
 	if(totalcost>oldtotalcost || (n>0 && totalcost==oldtotalcost)){
           fflush(NULL);
-	  fprintf(sp0,"Unexpected increase in total cost.  Breaking loop\n");
-	  break;
+	  fprintf(sp1,"Caution: Unexpected increase in total cost\n");
+	}
+        if(totalcost > mintotalcost){
+          nincreasedcostiter++;
+        }else{
+          nincreasedcostiter=0;
         }
       }
 
@@ -650,6 +663,12 @@ int UnwrapTile(infileT *infiles, outfileT *outfiles, paramT *params,
 
       /* find maximum flow on network, excluding arcs affected by masking */
       mostflow=MaxNonMaskFlow(flows,mag,nrow,ncol);
+      if(nincreasedcostiter>=mostflow){
+        fflush(NULL);
+        fprintf(sp0,"WARNING: Unexpected sustained increase in total cost."
+                "  Breaking loop\n");
+        break;
+      }
 
       /* break if we're done with all flow increments or problem is convex */
       if(nflowdone>=params->maxflow || nflowdone>=mostflow || params->p>=1.0){

src/snaphu.h

@@ -14,7 +14,7 @@
 /**********************/
 
 #define PROGRAMNAME          "snaphu"
-#define VERSION              "2.0.0"
+#define VERSION              "2.0.1"
 #define BUGREPORTEMAIL       "snaphu@gmail.com"
 #ifdef PI
 #undef PI
@@ -103,6 +103,7 @@
 #define DEF_VERBOSE          FALSE
 #define DEF_AMPLITUDE        TRUE
 #define AUTOCALCSTATMAX      0
+#define MAXNSHORTCYCLE       8192
 #define USEMAXCYCLEFRACTION  (-123)
 #define COMPLEX_DATA         1         /* file format */
 #define FLOAT_DATA           2         /* file format */

src/snaphu_cost.c

@@ -48,6 +48,13 @@ void **BuildStatCostsSmooth(float **wrappedphase, float **mag,
 			    long nrow, long ncol, tileparamT *tileparams, 
 			    outfileT *outfiles, paramT *params);
 static
+void MaskCost(costT *costptr);
+static
+void MaskSmoothCost(smoothcostT *smoothcostptr);
+static
+int MaskPrespecifiedArcCosts(void **costsptr, short **weights,
+                             long nrow, long ncol, paramT *params);
+static
 int GetIntensityAndCorrelation(float **mag, float **wrappedphase, 
                                float ***pwrptr, float ***corrptr, 
                                infileT *infiles, long linelen, long nlines,
@@ -200,12 +207,20 @@ int BuildCostArrays(void ***costsptr, short ***mstcostsptr,
 
   /* build or read the statistical cost arrays unless we were told not to */
   if(strlen(infiles->costinfile)){
+
+    /* read cost info from file */
     fprintf(sp1,"Reading cost information from file %s\n",infiles->costinfile);
     costs=NULL;
     Read2DRowColFile((void ***)&costs,infiles->costinfile,
 		     linelen,nlines,tileparams,costtypesize);
     (*costsptr)=costs;
 
+    /* weights of arcs next to masked pixels are set to zero */
+    /* make sure corresponding costs are nulled when costs are read from */
+    /*   file rather than internally generated since read costs are not */
+    /*   multiplied by weights */
+    MaskPrespecifiedArcCosts(costs,weights,nrow,ncol,params);
+
   }else if(params->costmode!=NOSTATCOSTS){
 
     /* get intensity and correlation info */
@@ -348,16 +363,23 @@ int BuildCostArrays(void ***costsptr, short ***mstcostsptr,
 	  tempcost=negcost;
 	}
 
-	/* clip scalar cost so it is between 0 and params->maxcost */
+	/* clip scalar cost so it is between 1 and params->maxcost */
+        /* note: weights used for MST algorithm will not be zero along */
+        /*   masked edges since they are clipped to 1, but MST is run */
+        /*   once on entire network, not just non-masked regions */
 	weights[row][col]=LClip(tempcost,MINSCALARCOST,params->maxcost);
 	
 	/* assign Lp costs if in Lp mode */
+        /* let scalar cost be zero if costs in both directions are zero */
 	if(params->p>=0){
 	  if(params->bidirlpn){
 	    bidircosts[row][col].posweight=LClip(poscost,0,params->maxcost);
 	    bidircosts[row][col].negweight=LClip(negcost,0,params->maxcost);
 	  }else{
 	    scalarcosts[row][col]=weights[row][col];
+            if(poscost==0 && negcost==0){
+              scalarcosts[row][col]=0;
+            }
 	  }
 	}
       }
@@ -582,10 +604,7 @@ void **BuildStatCostsTopo(float **wrappedphase, float **mag,
       if(colweight[row][col]==0){
 	  
 	/* masked pixel */
-	colcost[row][col].laycost=0;
-	colcost[row][col].offset=LARGESHORT/2;
-	colcost[row][col].dzmax=LARGESHORT;
-	colcost[row][col].sigsq=LARGESHORT;
+        MaskCost(&colcost[row][col]);
 
       }else{
 
@@ -732,10 +751,7 @@ void **BuildStatCostsTopo(float **wrappedphase, float **mag,
       if(rowweight[row][col]==0){
 	  
 	/* masked pixel */
-	rowcost[row][col].laycost=0;
-	rowcost[row][col].offset=LARGESHORT/2;
-	rowcost[row][col].dzmax=LARGESHORT;
-	rowcost[row][col].sigsq=LARGESHORT;
+        MaskCost(&rowcost[row][col]);
 
       }else{
 
@@ -898,10 +914,7 @@ void **BuildStatCostsDefo(float **wrappedphase, float **mag,
       if(colweight[row][col]==0){
 	  
 	/* masked pixel */
-	colcost[row][col].laycost=0;
-	colcost[row][col].offset=0;
-	colcost[row][col].dzmax=LARGESHORT;
-	colcost[row][col].sigsq=LARGESHORT;
+        MaskCost(&colcost[row][col]);
 
       }else{
 
@@ -964,10 +977,7 @@ void **BuildStatCostsDefo(float **wrappedphase, float **mag,
       if(rowweight[row][col]==0){
 	  
 	/* masked pixel */
-	rowcost[row][col].laycost=0;
-	rowcost[row][col].offset=0;
-	rowcost[row][col].dzmax=LARGESHORT;
-	rowcost[row][col].sigsq=LARGESHORT;
+        MaskCost(&rowcost[row][col]);
 
       }else{
 
@@ -1086,8 +1096,7 @@ void **BuildStatCostsSmooth(float **wrappedphase, float **mag,
       if(colweight[row][col]==0){
 	  
 	/* masked pixel */
-	colcost[row][col].offset=0;
-	colcost[row][col].sigsq=LARGESHORT;
+        MaskSmoothCost(&colcost[row][col]);
 
       }else{
 
@@ -1138,8 +1147,7 @@ void **BuildStatCostsSmooth(float **wrappedphase, float **mag,
       if(rowweight[row][col]==0){
 	  
 	/* masked pixel */
-	rowcost[row][col].offset=0;
-	rowcost[row][col].sigsq=LARGESHORT;
+        MaskSmoothCost(&rowcost[row][col]);
 
       }else{
 
@@ -1187,6 +1195,89 @@ void **BuildStatCostsSmooth(float **wrappedphase, float **mag,
 }
 
 
+/* function: MaskCost()
+ * --------------------
+ * Set values of costT structure pointed to by input pointer to give zero
+ * cost, as for arcs next to masked pixels.
+ */
+static
+void MaskCost(costT *costptr){
+
+  /* set to special values */
+  costptr->laycost=0;
+  costptr->offset=LARGESHORT/2;
+  costptr->dzmax=LARGESHORT;
+  costptr->sigsq=LARGESHORT;
+
+}
+
+
+/* function: MaskSmoothCost()
+ * --------------------------
+ * Set values of smoothcostT structure pointed to by input pointer to give zero
+ * cost, as for arcs next to masked pixels.
+ */
+static
+void MaskSmoothCost(smoothcostT *smoothcostptr){
+
+  /* set to special values */
+  smoothcostptr->offset=LARGESHORT/2;
+  smoothcostptr->sigsq=LARGESHORT;
+
+}
+
+
+/* function: MaskPrespecifiedArcCosts()
+ * ------------------------------------
+ * Loop over grid arcs and set costs to null if corresponding weights
+ * are null.
+ */
+static
+int MaskPrespecifiedArcCosts(void **costsptr, short **weights,
+                             long nrow, long ncol, paramT *params){
+
+  long row, col, maxcol;
+  costT **costs;
+  smoothcostT **smoothcosts;
+
+
+  /* set up pointers */
+  costs=NULL;
+  smoothcosts=NULL;
+  if(params->costmode==TOPO || params->costmode==DEFO){
+    costs=(costT **)costsptr;
+  }else if(params->costmode==SMOOTH){
+    smoothcosts=(smoothcostT **)costsptr;
+  }else{
+    fprintf(sp0,"illegal cost mode in MaskPrespecifiedArcCosts()\n");
+    exit(ABNORMAL_EXIT);
+  }
+
+  /* loop over all arcs */
+  for(row=0;row<2*nrow-1;row++){
+    if(row<nrow-1){
+      maxcol=ncol;
+    }else{
+      maxcol=ncol-1;
+    }
+    for(col=0;col<maxcol;col++){
+      if(weights[row][col]==0){
+        if(costs!=NULL){
+          MaskCost(&costs[row][col]);
+        }
+        if(smoothcosts!=NULL){
+          MaskSmoothCost(&smoothcosts[row][col]);
+        }
+      }
+    }
+  }
+
+  /* done */
+  return(0);
+  
+}
+
+
 /* function: GetIntensityAndCorrelation()
  * --------------------------------------
  * Reads amplitude and correlation info from files if specified.  If ampfile
@@ -1761,10 +1852,19 @@ void CalcCostTopo(void **costs, long flow, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((costT **)(costs))[arcrow][arccol];
-  dzmax=cost->dzmax;
   offset=cost->offset;
   sigsq=cost->sigsq;
+  dzmax=cost->dzmax;
   laycost=cost->laycost;
+
+  /* just return 0 if we have zero cost arc */
+  if(sigsq==LARGESHORT){
+    (*poscostptr)=0;
+    (*negcostptr)=0;
+    return;
+  }
+
+  /* compute argument to cost function */
   nshortcycle=params->nshortcycle;
   layfalloffconst=params->layfalloffconst;
   if(arcrow<nrow-1){
@@ -1864,6 +1964,15 @@ void CalcCostDefo(void **costs, long flow, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((costT **)(costs))[arcrow][arccol];
+
+  /* just return 0 if we have zero cost arc */
+  if(cost->sigsq==LARGESHORT){
+    (*poscostptr)=0;
+    (*negcostptr)=0;
+    return;
+  }
+
+  /* compute argument to cost function */
   nshortcycle=params->nshortcycle;
   layfalloffconst=params->layfalloffconst;
   idz1=labs(flow*nshortcycle+cost->offset);
@@ -1944,11 +2053,15 @@ void CalcCostSmooth(void **costs, long flow, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((smoothcostT **)(costs))[arcrow][arccol];
+
+  /* just return 0 if we have zero cost arc */
   if(cost->sigsq==LARGESHORT){
-    *poscostptr=0;
-    *negcostptr=0;
+    (*poscostptr)=0;
+    (*negcostptr)=0;
     return;
   }
+
+  /* compute argument to cost function */
   nshortcycle=params->nshortcycle;
   idz1=labs(flow*nshortcycle+cost->offset);
   idz2pos=labs((flow+nflow)*nshortcycle+cost->offset);
@@ -2238,6 +2351,24 @@ void CalcCostLPBiDir(void **costs, long flow, long arcrow, long arccol,
 /* function: CalcCostNonGrid()
  * ---------------------------
  * Calculates the arc cost given an array of long integer cost lookup tables.
+ *
+ * The cost array for each arc gives the cost for +/-flowmax units of
+ * flow around the flow value with minimum cost, which is not
+ * necessarily flow == 0.  The offset between the flow value with
+ * minimum cost and flow == 0 is given by arroffset = costarr[0].
+ * Positive flow values k for k = 1 to flowmax relative to this min
+ * cost flow value are in costarr[k].  Negative flow values k relative
+ * to the min cost flow from k = -1 to -flowmax costarr[flowmax-k].
+ * costarr[2*flowmax+1] contains a scaling factor for extrapolating
+ * beyond the ends of the cost table, assuming quadratically (with an offset)
+ * increasing cost (subject to rounding and scaling).
+ *
+ * As of summer 2019, the rationale for how seconeary costs are
+ * extrapolated beyond the end of the table has been lost to time, but
+ * the logic at least does give a self-consistent cost function that
+ * is continuous at +/-flowmax and quadratically increases beyond,
+ * albeit not necessarily with a starting slope that has an easily
+ * intuitive basis.
  */
 void CalcCostNonGrid(void **costs, long flow, long arcrow, long arccol, 
                      long nflow, long nrow, paramT *params, 
@@ -2345,6 +2476,13 @@ long EvalCostTopo(void **costs, short **flows, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((costT **)(costs))[arcrow][arccol];
+
+  /* just return 0 if we have zero cost arc */
+  if(cost->sigsq==LARGESHORT){
+    return(0);
+  }
+  
+  /* compute argument to cost function */
   if(arcrow<nrow-1){
 
     /* row cost: dz symmetric with respect to origin */
@@ -2388,6 +2526,13 @@ long EvalCostDefo(void **costs, short **flows, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((costT **)(costs))[arcrow][arccol];
+
+  /* just return 0 if we have zero cost arc */
+  if(cost->sigsq==LARGESHORT){
+    return(0);
+  }
+  
+  /* compute argument to cost function */
   idz1=labs(flows[arcrow][arccol]*(params->nshortcycle)+cost->offset);
 
   /* calculate and return cost */
@@ -2417,9 +2562,13 @@ long EvalCostSmooth(void **costs, short **flows, long arcrow, long arccol,
 
   /* get arc info */
   cost=&((smoothcostT **)(costs))[arcrow][arccol];
+
+  /* just return 0 if we have zero cost arc */
   if(cost->sigsq==LARGESHORT){
     return(0);
   }
+
+  /* compute argument to cost function */
   idz1=labs(flows[arcrow][arccol]*(params->nshortcycle)+cost->offset);
 
   /* calculate and return cost */

src/snaphu_io.c

@@ -899,6 +899,11 @@ int CheckParams(infileT *infiles, outfileT *outfiles,
     fprintf(sp0,"defomax exceeds range of short int for given nshortcycle\n");
     exit(ABNORMAL_EXIT);
   }
+  if(params->nshortcycle < 1 || params->nshortcycle > MAXNSHORTCYCLE){
+    fflush(NULL);
+    fprintf(sp0,"illegal value for nshortcycle\n");
+    exit(ABNORMAL_EXIT);
+  }
   if(params->maxnewnodeconst<=0 || params->maxnewnodeconst>1){
     fflush(NULL);
     fprintf(sp0,"maxnewnodeconst must be between 0 and 1\n");

src/snaphu_solver.c

@@ -2108,14 +2108,6 @@ void GetArcGrid(nodeT *from, nodeT *to, long *arcrow, long *arccol,
       *arccol=tocol;
       *arcdir=-1;
     }else{
-#define DIAG_GETARCGRID
-#ifdef DIAG_GETARCGRID
-      if(!(torow>0 && nodes[torow-1][tocol].group==BOUNDARYPTR)){
-        fflush(NULL);
-        fprintf(stderr,"BUG: should not have gotten here in GetArcGrid()\n");
-        exit(1);
-      }
-#endif
       *arcrow=torow+nrow-1;
       *arccol=tocol;
       *arcdir=1;
@@ -2134,14 +2126,6 @@ void GetArcGrid(nodeT *from, nodeT *to, long *arcrow, long *arccol,
       *arccol=fromcol;
       *arcdir=1;
     }else{
-#define DIAG_GETARCGRID
-#ifdef DIAG_GETARCGRID
-      if(!(fromrow>0 && nodes[fromrow-1][fromcol].group==BOUNDARYPTR)){
-        fflush(NULL);
-        fprintf(stderr,"BUG: should not have gotten here in GetArcGrid()\n");
-        exit(1);
-      }
-#endif
       *arcrow=fromrow+nrow-1;
       *arccol=fromcol;
       *arcdir=-1;

src/snaphu_tile.c

@@ -122,6 +122,8 @@ int SetRightEdge(long nrow, long ncol, long tilerow, long tilecol,
                  void **voidrightedgecosts, short **rightedgeflows, 
                  paramT *params, short **bulkoffsets);
 static
+short AvgSigSq(short sigsq1, short sigsq2);
+static
 int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes, 
                       nodesuppT **nodesupp, scndryarcT **scndryarcs, 
                       long ***scndrycosts, long *nnewnodesptr, 
@@ -1205,6 +1207,7 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
   long nrow, ncol, prevnrow, prevncol, nextnrow, nextncol;
   long n, ncycle, nflowdone, nflow, candidatelistsize, candidatebagsize;
   long nnodes, maxnflowcycles, arclen, narcs, sourcetilenum, flowmax;
+  long nincreasedcostiter;
   long *totarclens;
   long ***scndrycosts;
   double avgarclen;
@@ -1216,7 +1219,7 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
   short **tempregions, *regionsbelow, *regionsabove;
   int *nscndrynodes, *nscndryarcs;
   incrcostT **incrcosts;
-  totalcostT totalcost, oldtotalcost;
+  totalcostT totalcost, oldtotalcost, mintotalcost;
   nodeT *source;
   nodeT **scndrynodes, ***scndryapexes;
   signed char **iscandidate;
@@ -1370,8 +1373,8 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
 	}
 	scndrycosts[i][j][2*flowmax+1]=LRound(scndrycosts[i][j][2*flowmax+1]
 					      /avgarclen);
-	if(scndrycosts[i][j][2*flowmax+1]<1){
-	  scndrycosts[i][j][2*flowmax+1]=1;
+	if(scndrycosts[i][j][2*flowmax+1]<0){
+	  scndrycosts[i][j][2*flowmax+1]=0;
 	}
       }
     }
@@ -1414,7 +1417,9 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
 	      &candidatebagsize,&candidatebag,&candidatelistsize,
 	      &candidatelist,NULL,NULL,&bkts,&dummylong,NULL,NULL,NULL,
 	      NULL,NULL,NULL,NULL,ntiles,0,&notfirstloop,&totalcost,params);
-
+  oldtotalcost=totalcost;
+  mintotalcost=totalcost;
+  nincreasedcostiter=0;
 
   /* set pointers to functions for nongrid secondary network */
   CalcCost=CalcCostNonGrid;
@@ -1456,10 +1461,17 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
       oldtotalcost=totalcost;
       totalcost=EvaluateTotalCost((void **)scndrycosts,scndryflows,ntiles,0,
 				  nscndryarcs,params);
+      if(totalcost<mintotalcost){
+        mintotalcost=totalcost;
+      }
       if(totalcost>oldtotalcost || (n>0 && totalcost==oldtotalcost)){
         fflush(NULL);
-        fprintf(sp0,"Unexpected increase in total cost.  Breaking loop\n");
-        break;
+        fprintf(sp1,"Caution: Unexpected increase in total cost\n");
+      }
+      if(totalcost>mintotalcost){
+        nincreasedcostiter++;
+      }else{
+        nincreasedcostiter=0;
       }
     }
 
@@ -1471,6 +1483,15 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
       nflowdone=1;
     }
 
+    /* break if total cost increase is sustained */
+    if(nincreasedcostiter>=params->maxflow){
+      fflush(NULL);
+      fprintf(sp0,"WARNING: Unexpected sustained increase in total cost."
+              "  Breaking loop\n");
+      break;
+    }
+
+
     /* break if we're done with all flow increments or problem is convex */
     if(nflowdone>=params->maxflow){
       break;
@@ -2392,9 +2413,9 @@ int SetUpperEdge(long ncol, long tilerow, long tilecol, void **voidcosts,
 	dpsi-=1.0;
       }
       if(CalcCost==CalcCostTopo || CalcCost==CalcCostDefo){
-	upperedgecosts[0][col].offset=nshortcycle*dpsi;
-	upperedgecosts[0][col].sigsq=ceil((costs[0][col].sigsq
-					   +costsabove[col].sigsq)/2.0);
+	upperedgecosts[0][col].offset=(short )LRound(nshortcycle*dpsi);
+	upperedgecosts[0][col].sigsq=AvgSigSq(costs[0][col].sigsq,
+                                              costsabove[col].sigsq);
 	if(costs[0][col].dzmax>costsabove[col].dzmax){
 	  upperedgecosts[0][col].dzmax=costs[0][col].dzmax;
 	}else{
@@ -2406,9 +2427,9 @@ int SetUpperEdge(long ncol, long tilerow, long tilecol, void **voidcosts,
 	  upperedgecosts[0][col].laycost=costsabove[col].laycost;
 	}
       }else if(CalcCost==CalcCostSmooth){
-	upperedgesmoothcosts[0][col].offset=nshortcycle*dpsi;
-	upperedgesmoothcosts[0][col].sigsq=
-	  ceil((smoothcosts[0][col].sigsq+smoothcostsabove[col].sigsq)/2.0);
+	upperedgesmoothcosts[0][col].offset=(short )LRound(nshortcycle*dpsi);
+	upperedgesmoothcosts[0][col].sigsq
+          =AvgSigSq(smoothcosts[0][col].sigsq,smoothcostsabove[col].sigsq);
       }else if(CalcCost==CalcCostL0 || CalcCost==CalcCostL1
 	       || CalcCost==CalcCostL2 || CalcCost==CalcCostLP){
 	((short **)voidupperedgecosts)[0][col]=
@@ -2528,9 +2549,9 @@ int SetLowerEdge(long nrow, long ncol, long tilerow, long tilecol,
 	dpsi-=1.0;
       }
       if(CalcCost==CalcCostTopo || CalcCost==CalcCostDefo){
-	loweredgecosts[0][col].offset=nshortcycle*dpsi;
-	loweredgecosts[0][col].sigsq=ceil((costs[nrow-2][col].sigsq
-					   +costsbelow[col].sigsq)/2.0);
+	loweredgecosts[0][col].offset=(short )LRound(nshortcycle*dpsi);
+	loweredgecosts[0][col].sigsq=AvgSigSq(costs[nrow-2][col].sigsq,
+                                              costsbelow[col].sigsq);
 	if(costs[nrow-2][col].dzmax>costsbelow[col].dzmax){
 	  loweredgecosts[0][col].dzmax=costs[nrow-2][col].dzmax;
 	}else{
@@ -2542,10 +2563,9 @@ int SetLowerEdge(long nrow, long ncol, long tilerow, long tilecol,
 	  loweredgecosts[0][col].laycost=costsbelow[col].laycost;
 	}
       }else if(CalcCost==CalcCostSmooth){
-	loweredgesmoothcosts[0][col].offset=nshortcycle*dpsi;
-	loweredgesmoothcosts[0][col].sigsq=
-	  ceil((smoothcosts[nrow-2][col].sigsq
-		+smoothcostsbelow[col].sigsq)/2.0);
+	loweredgesmoothcosts[0][col].offset=(short )LRound(nshortcycle*dpsi);
+	loweredgesmoothcosts[0][col].sigsq
+	  =AvgSigSq(smoothcosts[nrow-2][col].sigsq,smoothcostsbelow[col].sigsq);
       }else if(CalcCost==CalcCostL0 || CalcCost==CalcCostL1
 	       || CalcCost==CalcCostL2 || CalcCost==CalcCostLP){
 	((short **)voidloweredgecosts)[0][col]=
@@ -2662,10 +2682,11 @@ int SetLeftEdge(long nrow, long prevncol, long tilerow, long tilecol,
 	dpsi-=1.0;
       }
       if(CalcCost==CalcCostTopo || CalcCost==CalcCostDefo){
-	leftedgecosts[row][0].offset=(TILEDPSICOLFACTOR*nshortcycle*dpsi);
-	leftedgecosts[row][0].sigsq=
-	  ceil((costs[row+nrow-1][0].sigsq
-		+lastcosts[row+nrow-1][prevncol-2].sigsq)/2.0);
+	leftedgecosts[row][0].offset=(short )LRound(TILEDPSICOLFACTOR
+                                                    *nshortcycle*dpsi);
+	leftedgecosts[row][0].sigsq
+          =AvgSigSq(costs[row+nrow-1][0].sigsq,
+                    lastcosts[row+nrow-1][prevncol-2].sigsq);
 	if(costs[row+nrow-1][0].dzmax>lastcosts[row+nrow-1][prevncol-2].dzmax){
 	  leftedgecosts[row][0].dzmax=costs[row+nrow-1][0].dzmax;
 	}else{
@@ -2680,10 +2701,10 @@ int SetLeftEdge(long nrow, long prevncol, long tilerow, long tilecol,
 	}
       }else if(CalcCost==CalcCostSmooth){
 	leftedgesmoothcosts[row][0].offset
-	  =(TILEDPSICOLFACTOR*nshortcycle*dpsi);
-	leftedgesmoothcosts[row][0].sigsq=
-	  ceil((smoothcosts[row+nrow-1][0].sigsq
-		+lastsmoothcosts[row+nrow-1][prevncol-2].sigsq)/2.0);
+	  =(short )LRound(TILEDPSICOLFACTOR*nshortcycle*dpsi);
+	leftedgesmoothcosts[row][0].sigsq
+	  =AvgSigSq(smoothcosts[row+nrow-1][0].sigsq,
+                    lastsmoothcosts[row+nrow-1][prevncol-2].sigsq);
       }else if(CalcCost==CalcCostL0 || CalcCost==CalcCostL1
 	       || CalcCost==CalcCostL2 || CalcCost==CalcCostLP){
 	((short **)voidleftedgecosts)[row][0]=
@@ -2807,10 +2828,11 @@ int SetRightEdge(long nrow, long ncol, long tilerow, long tilecol,
 	dpsi-=1.0;
       }
       if(CalcCost==CalcCostTopo || CalcCost==CalcCostDefo){
-	rightedgecosts[row][0].offset=(TILEDPSICOLFACTOR*nshortcycle*dpsi);
+	rightedgecosts[row][0].offset=(short )LRound(TILEDPSICOLFACTOR
+                                                     *nshortcycle*dpsi);
 	rightedgecosts[row][0].sigsq
-	  =ceil((costs[row+nrow-1][ncol-2].sigsq
-		 +nextcosts[row+nrow-1][0].sigsq)/2.0);
+	  =AvgSigSq(costs[row+nrow-1][ncol-2].sigsq,
+                    nextcosts[row+nrow-1][0].sigsq);
 	if(costs[row+nrow-1][ncol-2].dzmax>nextcosts[row+nrow-1][0].dzmax){
 	  rightedgecosts[row][0].dzmax=costs[row+nrow-1][ncol-2].dzmax;
 	}else{
@@ -2823,10 +2845,10 @@ int SetRightEdge(long nrow, long ncol, long tilerow, long tilecol,
 	}
       }else if(CalcCost==CalcCostSmooth){
 	rightedgesmoothcosts[row][0].offset
-	  =(TILEDPSICOLFACTOR*nshortcycle*dpsi);
+	  =(short )LRound(TILEDPSICOLFACTOR*nshortcycle*dpsi);
 	rightedgesmoothcosts[row][0].sigsq
-	  =ceil((smoothcosts[row+nrow-1][ncol-2].sigsq
-		 +nextsmoothcosts[row+nrow-1][0].sigsq)/2.0);
+	  =AvgSigSq(smoothcosts[row+nrow-1][ncol-2].sigsq,
+                    nextsmoothcosts[row+nrow-1][0].sigsq);
       }else if(CalcCost==CalcCostL0 || CalcCost==CalcCostL1
 	       || CalcCost==CalcCostL2 || CalcCost==CalcCostLP){
 	((short **)voidrightedgecosts)[row][0]=
@@ -2908,6 +2930,34 @@ int SetRightEdge(long nrow, long ncol, long tilerow, long tilecol,
 }
 
 
+/* function: AvgSigSq()
+ * --------------------
+ * Return average of sigsq values after chcking for special value and
+ * clipping to short.
+ */
+static
+short AvgSigSq(short sigsq1, short sigsq2){
+
+  int sigsqavg;
+
+  
+  /* if either value is special LARGESHORT value, use that */
+  if(sigsq1==LARGESHORT || sigsq2==LARGESHORT){
+    return(LARGESHORT);
+  }
+
+  /* compute average */
+  sigsqavg=(int )ceil(0.5*(((int )sigsq1)+((int )sigsq2)));
+
+  /* clip */
+  sigsqavg=LClip(sigsqavg,-LARGESHORT,LARGESHORT);
+
+  /* return */
+  return((short )sigsqavg);
+  
+}
+
+
 /* function: TraceSecondaryArc()
  * -----------------------------
  */
@@ -2931,7 +2981,7 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
   long i, row, col, nnewnodes, arclen, ntilerow, ntilecol, arcnum;
   long tilenum, nflow, primaryarcrow, primaryarccol, poscost, negcost, nomcost;
   long nnrow, nncol, calccostnrow, nnewarcs, arroffset, nshortcycle;
-  long mincost, mincostflow, minweight;
+  long mincost, mincostflow, minweight, maxcost;
   long *scndrycostarr;
   double sigsq, sumsigsqinv, tempdouble, tileedgearcweight;
   short **flows;
@@ -3096,6 +3146,9 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
 
       /* keep absolute cost of arc to the previous node */
       if(!zerocost){
+
+        /* accumulate incremental cost in table for each nflow increment */
+        /* offset flow in flow array temporarily by arroffset then undo below */
 	flows[primaryarcrow][primaryarccol]-=primaryarcdir*arroffset;
 	nomcost=EvalCost(costs,flows,primaryarcrow,primaryarccol,calccostnrow,
 			 params);
@@ -3125,23 +3178,36 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
 	  }
 	}
 	flows[primaryarcrow][primaryarccol]+=primaryarcdir*arroffset;
+
+        /* accumulate term to be used for cost growth beyond table bounds */
 	if(CalcCost==CalcCostTopo || CalcCost==CalcCostDefo){
 	  sigsq=((costT **)costs)[primaryarcrow][primaryarccol].sigsq;
 	}else if(CalcCost==CalcCostSmooth){
 	  sigsq=((smoothcostT **)costs)[primaryarcrow][primaryarccol].sigsq;
 	}else if(CalcCost==CalcCostL0 || CalcCost==CalcCostL1 
 		 || CalcCost==CalcCostL2 || CalcCost==CalcCostLP){
-	  sigsq=((short **)costs)[primaryarcrow][primaryarccol];
+	  minweight=((short **)costs)[primaryarcrow][primaryarccol];
+          if(minweight<1){
+            sigsq=LARGESHORT;
+          }else{
+            sigsq=1.0/(double )minweight;
+          }
 	}else if(CalcCost==CalcCostL0BiDir || CalcCost==CalcCostL1BiDir
 		 || CalcCost==CalcCostL2BiDir || CalcCost==CalcCostLPBiDir){
 	  minweight=LMin(((bidircostT **)costs)[primaryarcrow][primaryarccol]
 			 .posweight,
 			 ((bidircostT **)costs)[primaryarcrow][primaryarccol]
 			 .negweight);
+          if(minweight<1){
+            sigsq=LARGESHORT;
+          }else{
             sigsq=1.0/(double )minweight;
           }
+	}
+        if(sigsq<LARGESHORT){    /* zero cost arc if sigsq == LARGESHORT */
           sumsigsqinv+=(1.0/sigsq);
         }
+      }
 
       /* break if found the secondary arc tail */
       if(primarytail->group==ONTREE){
@@ -3161,6 +3227,7 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
 
     /* find flow index with minimum cost */
     mincost=0;
+    maxcost=0;
     mincostflow=0;
     for(nflow=1;nflow<=flowmax;nflow++){
       if(scndrycostarr[nflow]<mincost){
@@ -3171,6 +3238,18 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
 	mincost=scndrycostarr[flowmax+nflow];
 	mincostflow=-nflow;
       }
+      if(scndrycostarr[nflow]>maxcost){
+	maxcost=scndrycostarr[nflow];
+      }
+      if(scndrycostarr[flowmax+nflow]>maxcost){
+	maxcost=scndrycostarr[flowmax+nflow];
+      }
+    }
+
+    /* if cost was all zero, treat as zero cost arc */
+    if(maxcost==mincost){
+      zerocost=TRUE;
+      sumsigsqinv=0;
     }
 
     /* break if cost array adequately centered on minimum cost flow */
@@ -3198,12 +3277,16 @@ int TraceSecondaryArc(nodeT *primaryhead, nodeT **scndrynodes,
     return(0);
   }
 
-
   /* see if we have a secondary arc on the edge of the full-sized array */
   /* these arcs have zero cost since the edge is treated as a single node */
+  /* secondary arcs whose primary arcs all have zero cost are also zeroed */
   if(zerocost){
 
     /* set sum of standard deviations to indicate zero-cost secondary arc */
+    scndrycostarr[0]=0;
+    for(nflow=1;nflow<=2*flowmax;nflow++){
+      scndrycostarr[nflow]=0;
+    }
     scndrycostarr[2*flowmax+1]=ZEROCOSTARC;
 
   }else{
@@ -4042,12 +4125,6 @@ int AssembleTileConnComps(long linelen, long nlines,
               conncompsizes[nconncomp].icompfull=0;
               conncompsizes[nconncomp].npix=tileconncompsizes[k].npix;
               nconncomp++;
-#define DEBUG              
-#ifdef DEBUG              
-if(nconncomp>nconncompmem){
-  fprintf(sp0,"ERROR--THIS IS A BUG\n");
-}
-#endif
             }
           }