Correct bin size for histogram of radiometric indices (-1.0 to 1.0 valued)

src/rasterstats1d.cpp

@@ -159,10 +159,12 @@ void RasterStats1D::computeMean(int rn) {
 	sums = new float*[objMap->nObjects];
 	means = new float*[objMap->nObjects];
 	for (t=0;t<objMap->nObjects;t++) {
-		counts[t] = new float[1];
-		counts[t][0] = 0.0;
+		counts[t] = new float[B];
 		sums[t] = new float[B];
-		for (k=0;k<B;k++) sums[t][k] = 0.0;
+		for (k=0;k<B;k++) {
+			sums[t][k] = 0.0;
+			counts[t][k] = 0.0;
+		}
 		means[t] = new float[B];
 	}
 
@@ -179,17 +181,19 @@ void RasterStats1D::computeMean(int rn) {
 					if (d != rasters[rn]->getNoValue(k)) {
 						objIdx = objMap->rLabels[curr];
 						sums[objIdx][k] += w*d;
-						if (k==0) counts[objIdx][0] += w;
+						counts[objIdx][k] += w;
 					}
 				}
 			}
 		}
 		for (t=0;t<objMap->nObjects;t++) {
-			means[t][k] = sums[t][k] / counts[t][0];
+			if (counts[t][k] > 0)
+				means[t][k] = sums[t][k] / counts[t][k];
+			else means[t][k] = rasters[rn]->getNoValue(k);
 		}
 	}
 
-	stats_sizes[rn][OSTAT_COUNT] = 1;
+	stats_sizes[rn][OSTAT_COUNT] = B;
 	stats[rn][OSTAT_COUNT] = counts;
 	stats_sizes[rn][OSTAT_SUM] = B;
 	stats[rn][OSTAT_SUM] = sums;
@@ -233,7 +237,9 @@ void RasterStats1D::computeStd(int rn) {
 			}
 		}
 		for (t=0;t<objMap->nObjects;t++) {
-			stds[t][k] = sqrt(sums[t][k] / stats[rn][OSTAT_COUNT][t][0]);
+			if (stats[rn][OSTAT_COUNT][t][k] > 0)
+				stds[t][k] = sqrt(sums[t][k] / stats[rn][OSTAT_COUNT][t][k]);
+			else stds[t][k] = rasters[rn]->getNoValue(k);
 		}
 	}
 
@@ -356,17 +362,21 @@ void RasterStats1D::computeHistograms(int rn, float bs) {
 	Data ***bounds = new Data**[objMap->nObjects];
 	long ***tmp = new long**[objMap->nObjects];
 
+	float bsr;
+
 	for (t=0;t<objMap->nObjects;t++) {
 		tmp[t] = new long*[B];
 		bounds[t] = new Data*[B];
 		R[t] = new int[B];
 		for (k=0;k<B;k++) {
-			R[t][k] = (int)round((float)(stats[rn][OSTAT_MAX][t][k] - stats[rn][OSTAT_MIN][t][k] + 1) / bs);
+			// Workaround for non-scaled radiometric indices (-1.0 to 1.0 valued)
+			bsr = bs < ((stats[rn][OSTAT_MAX][t][k] - stats[rn][OSTAT_MIN][t][k]) / 100) ? bs : bs/1000;
+			R[t][k] = (int)round((float)(stats[rn][OSTAT_MAX][t][k] - stats[rn][OSTAT_MIN][t][k] + bsr) / bsr);
 			tmp[t][k] = new long[R[t][k]+1]; // idx R for total
 			bounds[t][k] = new Data[R[t][k]+1]; // R bins --> R+1 bounds
-			minH = (float)(stats[rn][OSTAT_MIN][t][k])-bs/2;
+			minH = (float)(stats[rn][OSTAT_MIN][t][k])-bsr/2;
 			for (s=0;s<=R[t][k];s++) {
-				bounds[t][k][s] = minH + s * bs;
+				bounds[t][k][s] = minH + s * bsr;
 				tmp[t][k][s] = 0;
 			}
 		}