Plot text values for partial pressure graphs

The algorithms attempt to identify "interesting" points where the user
might want to know the value of the graph.

Signed-off-by: Dirk Hohndel <dirk@hohndel.org>

profile.c

@@ -8,6 +8,7 @@
 #include <stdarg.h>
 #include <string.h>
 #include <time.h>
+#include <math.h>
 
 #include "dive.h"
 #include "display.h"
@@ -505,6 +506,264 @@ static void plot_depth_scale(struct graphics_context *gc, struct plot_info *pi)
 	}
 }
 
+/* ap points to an array of int with pi->nr + 1 elements that is
+ * ininitialized with just one -1 entry
+ * this adds entries (if they aren't too close to an existing one)
+ * and keeps things sorted
+ * we KNOW the array is big enough to hold all possible indices
+ * a2p is a secondary array - we insert value at the same relative
+ * positio as idx in ap */
+static void add_index(int idx, int margin, int **ap, int **a2p, int value)
+{
+	int j, i = 0;
+	int *a = *ap;
+	int *a2 = *a2p;
+
+	while (a[i] != -1 && a[i] < idx)
+		i++;
+	if (a[i] == idx)
+		return;
+	if (a[i] != -1 && a[i - 1] != -1 && idx - a[i - 1] < margin)
+		return;
+	if (a[i] != -1 && a[i] - idx < margin)
+		return;
+	j = i;
+	while (a[j] != -1)
+		j++;
+	while (j >= i) {
+		a[j+1] = a[j];
+		a2[j+1] = a2[j];
+		j--;
+	}
+	a[i] = idx;
+	a2[i] = value;
+}
+
+#define LI(_i,_j) MAX((_i)-(_j), 0)
+#define RI(_i,_j) MIN((_i)+(_j), nr - 1)
+#define SPIKE(_i,_s) if (fabs(_s) > fabs(spk_data[_i])) spk_data[_i] = (_s)
+/* this is an attempt at a metric that finds spikes in a data series */
+static void calculate_spikyness(int nr, double *data, double *spk_data, int deltax, double deltay)
+{
+	int i, j;
+	double dminl, dminr, dmaxl, dmaxr;
+
+#if DEBUG_PROFILE > 2
+	printf("Spike data: \n 0 ");
+#endif
+	for (i = 0; i < nr; i++) {
+		dminl = dminr = dmaxl = dmaxr = data[i];
+		spk_data[i] = 0.0;
+		for (j = 1; j < deltax; j++) {
+			if (data[LI(i,j)] < dminl)
+				dminl = data[LI(i,j)];
+			if (data[LI(i,j)] > dmaxl)
+				dmaxl = data[LI(i,j)];
+
+			if (data[RI(i,j)] < dminr)
+				dminr = data[RI(i,j)];
+			if (data[RI(i,j)] > dmaxr)
+				dmaxr = data[RI(i,j)];
+
+			/* don't do super narrow */
+			if (j < deltax / 3)
+				continue;
+			/* falling edge on left */
+			if (dmaxl == data[i] && dmaxr - data[i] < 0.1 * (data[i] - dminl))
+				SPIKE(i,(data[i] - dminl) / j);
+			/* falling edge on right */
+			if (dmaxr == data[i] && dmaxl - data[i] < 0.1 * (data[i] - dminr))
+				SPIKE(i,(data[i] - dminr) / j);
+
+			/* minima get a negative spike value */
+			/* rising edge on left */
+			if (dminl == data[i] && data[i] - dminr < 0.1 * (dmaxl - data[i]))
+				SPIKE(i,(data[i] - dmaxl) / j);
+			/* rising edge on right */
+			if (dminr == data[i] && data[i] - dminl < 0.1 * (dmaxr - data[i]))
+				SPIKE(i,(data[i] - dmaxr) / j);
+		}
+#if DEBUG_PROFILE > 2
+		fprintf(debugfile, "%.4lf ", spk_data[i]);
+		if (i % 12 == 11)
+			fprintf(debugfile, "\n%2d ", (i + 1) / 12);
+#endif
+	}
+#if DEBUG_PROFILE > 2
+	printf("\n");
+#endif
+}
+
+/* only show one spike in a deltax wide region - pick the highest (and first if the same) */
+static gboolean higher_spike(double *spk_data, int idx, int nr, int deltax)
+{
+	int i;
+	double s = fabs(spk_data[idx]);
+	for (i = MAX(0, idx - deltax); i <= MIN(idx + deltax, nr - 1); i++)
+		if (fabs(spk_data[i]) > s)
+			return TRUE;
+		else if (fabs(spk_data[i]) == s && i < idx)
+			return TRUE;
+	return FALSE;
+}
+
+/* this figures out which time stamps provide "interesting" formations in the graphs;
+ * this includes local minima and maxima as well as long plateaus.
+ * pass in the function that returns the value at a certain point (as double),
+ * the delta in time (expressed as number of data points of "significant time")
+ * the delta at which the value is considered to have been "significantly changed" and
+ * the number of points to cover
+ * returns a list of indices that ends with a -1 of times that are "interesting" */
+static void find_points_of_interest(struct plot_info *pi, double (*value_func)(int, struct plot_info *),
+				int deltax, double deltay, int **poip, int **poip_vpos)
+{
+	int i, j, nr = pi->nr;
+	double *data, *data_max, *data_min, *spk_data;
+	double min, max;
+	int *pois;
+
+	/* avoid all the function calls by creating a local array and
+	 * have some helper arrays to make our lifes easier */
+
+	data = malloc(nr * sizeof(double));
+	data_max = malloc(nr * sizeof(double));
+	data_min = malloc(nr * sizeof(double));
+	spk_data = malloc(nr * sizeof(double));
+
+	pois = *poip = malloc((nr + 1) * sizeof(int));
+	*poip_vpos = malloc((nr + 1) * sizeof(int));
+	pois[0] = -1;
+	pois[1] = -1;
+
+	/* copy the data and get the absolute minimum and maximum while we do it */
+	for (i = 0; i < nr; i++) {
+		data_max[i] = data_min[i] = data[i] = value_func(i, pi);
+		if (i == 0 || data[i] < min)
+			min = data[i];
+		if (i == 0 || data[i] > max)
+			max = data[i];
+	}
+	/* next find out if there are real spikes in the graph */
+	calculate_spikyness(nr, data, spk_data, deltax, deltay);
+
+	/* now process all data points */
+	for (i = 0; i < nr; i++) {
+		/* get the local min/max */
+		for (j = MAX(0, i - deltax); j < i + deltax && j < nr; j++) {
+			if (data[j] < data[i])
+				data_min[i] = data[j];
+			if (data[j] > data[i])
+				data_max[i] = data[j];
+		}
+		/* is i the overall minimum or maximum */
+		if (data[i] == max)
+			add_index(i, deltax, poip, poip_vpos, BOTTOM);
+		if (data[i] == min)
+			add_index(i, deltax, poip, poip_vpos, TOP);
+		/* is i a spike? */
+		if (fabs(spk_data[i]) > 0.01 && ! higher_spike(spk_data, i, nr, deltax)) {
+			if (spk_data[i] > 0.0)
+				add_index(i, deltax, poip, poip_vpos, BOTTOM);
+			if (spk_data[i] < 0.0)
+				add_index(i, deltax, poip, poip_vpos, TOP);
+		}
+		/* is i a significant local minimum or maximum? */
+		if (data[i] == data_min[i] && data_max[i] - data[i] > deltay)
+			add_index(i, deltax, poip, poip_vpos, TOP);
+		if (data[i] == data_max[i] && data[i] - data_min[i] > deltay)
+			add_index(i, deltax, poip, poip_vpos, BOTTOM);
+	}
+	/* still need to search for plateaus */
+}
+
+static void setup_pp_limits(struct graphics_context *gc, struct plot_info *pi)
+{
+	int maxdepth;
+
+	gc->leftx = 0;
+	gc->rightx = get_maxtime(pi);
+
+	/* the maxdepth already includes extra vertical space - and if
+	 * we take the corresponding pressure as maximum partial
+	 * pressure the graph seems to look fine*/
+	maxdepth = get_maxdepth(pi);
+	gc->topy = (maxdepth + 10000) / 10000.0 * 1.01325;
+	gc->bottomy = 0.0;
+}
+
+static void plot_single_pp_text(struct graphics_context *gc, int sec, double pp, color_indice_t color)
+{
+	text_render_options_t tro = {12, color, CENTER, BOTTOM};
+	plot_text(gc, &tro, sec, pp, "%.1lf", pp);
+
+	if (color == PN2)
+		printf("pN2 %lf\n", pp);
+}
+
+#define MAXPP(_mpp, _pp) { _mpp = 0;			\
+	for(i = 0; i< pi->nr; i++)			\
+		if (pi->entry[i]._pp > _mpp)		\
+			_mpp = pi->entry[i]._pp;	\
+	}
+
+static double po2_value(int idx, struct plot_info *pi)
+{
+	return pi->entry[idx].po2;
+}
+
+static double pn2_value(int idx, struct plot_info *pi)
+{
+	return pi->entry[idx].pn2;
+}
+
+static double phe_value(int idx, struct plot_info *pi)
+{
+	return pi->entry[idx].phe;
+}
+
+static void plot_single_gas_pp_text(struct graphics_context *gc, struct plot_info *pi,
+				double (*value_func)(int, struct plot_info *),
+				double value_threshold, int color)
+{
+	int *pois, *pois_vpos;
+	int i, two_minutes = 1;
+	/* don't bother with local min/max if the dive is under two minutes */
+	if (pi->entry[pi->nr - 1].sec > 120) {
+		int idx = 0;
+		while (pi->entry[idx].sec == 0)
+			idx++;
+		while (pi->entry[idx + two_minutes].sec < 120)
+			two_minutes++;
+	} else {
+		two_minutes = pi->nr;
+	}
+	find_points_of_interest(pi, value_func, two_minutes, value_threshold, &pois, &pois_vpos);
+	for (i = 0; pois[i] != -1; i++) {
+		struct plot_data *entry = pi->entry + pois[i];
+#if DEBUG_PROFILE > 1
+		fprintf(debugfile, "POI at %d sec value %lf\n", entry->sec, entry->po2);
+#endif
+		plot_single_pp_text(gc, entry->sec, value_func(pois[i], pi), color);
+	}
+	free(pois);
+	free(pois_vpos);
+}
+
+static void plot_pp_text(struct graphics_context *gc, struct plot_info *pi)
+{
+	setup_pp_limits(gc, pi);
+
+	if (enabled_graphs.po2) {
+		plot_single_gas_pp_text(gc, pi, po2_value, 0.4, PO2);
+	}
+	if (enabled_graphs.pn2) {
+		plot_single_gas_pp_text(gc, pi, pn2_value, 0.4, PN2);
+	}
+	if (enabled_graphs.phe) {
+		plot_single_gas_pp_text(gc, pi, phe_value, 0.4, PHE);
+	}
+}
+
 static void plot_pp_gas_profile(struct graphics_context *gc, struct plot_info *pi)
 {
 	int i;
@@ -1591,8 +1850,10 @@ void plot(struct graphics_context *gc, cairo_rectangle_t *drawing_area, struct d
 	cairo_close_path(gc->cr);
 	cairo_stroke(gc->cr);
 
-	if (GRAPHS_ENABLED)
+	if (GRAPHS_ENABLED) {
 		plot_pp_gas_profile(gc, pi);
+		plot_pp_text(gc, pi);
+	}
 
 	/* now shift the translation back by half the margin;
 	 * this way we can draw the vertical scales on both sides */