[Git][debian-gis-team/snaphu][upstream] New upstream version 2.0.1
Antonio Valentino
gitlab at salsa.debian.org
Sun Sep 8 12:11:45 BST 2019
Antonio Valentino pushed to branch upstream at Debian GIS Project / snaphu
Commits:
c7fff519 by Antonio Valentino at 2019-09-08T11:00:06Z
New upstream version 2.0.1
- - - - -
6 changed files:
- src/snaphu.c
- src/snaphu.h
- src/snaphu_cost.c
- src/snaphu_io.c
- src/snaphu_solver.c
- src/snaphu_tile.c
Changes:
=====================================
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,¬firstloop,&totalcost,params);
+ oldtotalcost=totalcost;
+ mintotalcost=totalcost;
+ nincreasedcostiter=0;
/* regrow regions with -G parameter */
if(params->regrowconncomps){
@@ -633,11 +639,18 @@ 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;
+ }
}
/* consider this flow increment done if not too many neg cycles found */
@@ -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 at 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;
@@ -1227,7 +1230,7 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
bucketT *bkts;
char filename[MAXSTRLEN];
-
+
/* set up */
fprintf(sp1,"Assembling tiles\n");
ntilerow=params->ntilerow;
@@ -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;
}
}
}
@@ -1408,13 +1411,15 @@ int AssembleTiles(outfileT *outfiles, paramT *params,
incrcosts[i]=(incrcostT *)MAlloc(nscndryarcs[i]*sizeof(incrcostT));
nnodes+=nscndrynodes[i];
}
-
+
/* set up network for secondary solver */
InitNetwork(scndryflows,&dummylong,&ncycle,&nflowdone,&dummylong,&nflow,
&candidatebagsize,&candidatebag,&candidatelistsize,
&candidatelist,NULL,NULL,&bkts,&dummylong,NULL,NULL,NULL,
NULL,NULL,NULL,NULL,ntiles,0,¬firstloop,&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,22 +3178,35 @@ 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);
- sigsq=1.0/(double )minweight;
+ if(minweight<1){
+ sigsq=LARGESHORT;
+ }else{
+ sigsq=1.0/(double )minweight;
+ }
}
- sumsigsqinv+=(1.0/sigsq);
+ if(sigsq<LARGESHORT){ /* zero cost arc if sigsq == LARGESHORT */
+ sumsigsqinv+=(1.0/sigsq);
+ }
}
/* break if found the secondary arc tail */
@@ -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
}
}
View it on GitLab: https://salsa.debian.org/debian-gis-team/snaphu/commit/c7fff519947050f80cd102a77b25912796a1099d
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/snaphu/commit/c7fff519947050f80cd102a77b25912796a1099d
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20190908/9f2f3e9c/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list