Using polynomial interpolation to extract numerical data from seismic monitor images

Translated and modified from original article on Qiita. Disclaimer: I am not the author of the original article.

Hi everyone, I'm Wolf and I'm an Android, backend, web and desktop developer. This is my first article on Dev.to, so I hope you enjoy. I decided to translate the above linked article because I find it interesting to me, and it is also helpful for a project I'm working on. I've also slightly modified this article, so there will be additional code snippets and things which aren't present in the original article. The "seismic data" and all seismics, earthquakes related things here are for Japan.

This article will be a bit long, but I hope it isn't too confusing. (you can jump to the code section near the end of the article if you don't want to read the lengthy explanations)

Introduction

In third-party software that uses seismic monitor data, a specific number of features rely on the utilization of numerical values derived from GIF images. In fact, since the information is encoded in specific colors at specific pixel locations, the software needs a way to convert arbitrary colors into numeric values as well as a database to associate pixels with seismic observation points.

Applicable uses are the display of measured seismic intensity values or earthquake detection algorithms. Several types of data can be obtained, including real-time seismic intensity, maximum ground acceleration (PGA), maximum ground velocity, each using a different scale.

At the time of writing, the Japanese National Institute for Earth Science and Disaster Prevention (NIED) has not disclosed any official way to find the numerical values or a REST/WebSocket API (e.g. in JSON format) that features the numerical values of the data channels, so we need to rely on other methods.

There are several ways to get the numerical values, such as using tables that maps numerical data with red, green and blue (R, G, B) values. This article describes how to use polynomial interpolation to convert any color to any value on the scale. This is used in seismic monitor software to create navigable maps and to utilize numerical values for various functions.

Popular seismic monitor software that uses this method include JQuake and Kyoshin EEW Viewer for ingen.

Goal

During the development of the EEW Event Viewer software (old name of JQuake software, see blog about JQuake development history in Japanese here), a way to extract data from the following image formats is needed:

• GIF for real-time data
• PNG for old data from the server (at that time old data was still available, e.g. for the 2011 Great East Japan Earthquake)
• AVI videos of the data converted to PNG frames from the NIED data website

Considering the fact that colors vary slightly from format to format due to compression, I chose to use polynomial interpolation instead of using color mapping tables to extract numerical data.

What is polynomial interpolation?

If the function expression is unknown but the input and output is known, one way to approximate the function expression is to use polynomial interpolation. The goal is to be able to provide an output for future inputs to the unknown expression.

For example, in the following figure, there are seven known outputs from integer inputs (red) to a particular function. However, we want to know the output value 3.55 (green). This is the case when we want to get slightly off color values (for example, because of compression algorithms).

Performing polynomial interpolation on the red points yields the coefficients {an,...,a0} of the optimal polynomial function, which is expressed as:

where n is the degree of the polynomial. Here, a polynomial of degree 7 is used, and its curve is shown in blue.

Now we can know the output of 3.55 by entering 3.55 into the polynomial we found. The result is -0.4.

The example presented was a case of perfect agreement between a polynomial and a small set of points. However, if much more points are met (>50), using a very high polynomial becomes resource intensive, so usually keep the degree below 10 and accept the tolerances. Splitting the input into separate intervals when aspects of the unknown function change significantly is also a good solution.

Formulation of the problem

Looking back at the problem, we can see that all the data received on any channel has one point in common. In other words, the color scale is always the same. Therefore, our goal is to know a function that converts a color (R, G, B) to a position on the scale (say from 0 to 1). From this position, we can compute any value for the desired data channel, knowing the equation that converts the position on the scale to the type of data received.

(For dark theme users: If you have trouble seeing the image, please switch to light theme)

Polynomial interpolation

As mentioned before, we need to find a function that converts colors to positions. For this purpose, an analysis of the color scale image found in the seismic monitor website is done. Matlab was used for analysis. Vertical lines spanning from pixel (6,18) to (6,296) are extracted. Then there are 279 different colors for interpolation, pixel (6,296) corresponds to position 0.0 and (6,18) corresponds to position 1.0.

Using R, G, B

Plotting the position as a function of red, green and blue components of the extracted pixel lines yields the following figure:

To fit a polynomial function, it's important that one value of the abscissa gives at most one value of the ordinate. Otherwise, it cannot be called a function. We can see that the green component is not a function, while the red function is clearly defined up to about 180 values (covering the position from 0.22 to 0.4). For the blue component (covering the position from 0.15 to 0.4) it is up to 210. Also note that the outlier around 180 is due to the scale. In conclusion, it is not possible to perform polynomial interpolation on the R, G, and B values to get the positions.

Using H, S and V

You may have noticed that the color in the scale evolves in a linear fashion (from blue to red), just like the colors of the rainbow. In fact, it's evolving in strictly decreasing hues. The secret is to convert RGB to HSV color space to get an exploitable function.

If it is greater than 1, we can see that the H component (hue) is clearly defined and covers position values from 0 to 0.9. Fortunately, when the H component is less than 1, the V component is clearly defined and covers position values from 0.9 to 1.

Note that the rest of the article will consider the values of H, S and V to be between 0 and 1, not 0 to 255. This is explained by the fact that Java and Matlab return such values. Each of the following interpolations is performed with these values in mind. There will be a C++ snippet later in addition to the Java snippet in the original article, for people using languages other than Java.

First, the below are the results of the interpolation using the part of the H component covering the position from 0 to 0.6:

On the left side of the figure are computed coefficients (p_1 to p_7, corresponding to a6 to a0 in the previous equation). The subfigure on the upper right shows the set of points with the approximated polynomials. A polynomial of degree 6 has been selected and the red outliers have been excluded. We can see that the residuals (bottom right subfigure) are contained between -0.0075 and 0.008.

Next, the part of the H component covering the position from 0.6 to 0.9 is interpolated:

It turns out that the polynomial of degree 4 is sufficient to have a residual error of less than 0.004 in absolute value.

Finally, the V values covering the positions between 0.9 and 1.0 are covered in the third interpolation:

When outliers are excluded and the degree of the polynomial is set to 2, the residuals become very small (less than 0.001).

Subsequent work has shown discrepancies with values measured in historical seismic data at the NIED site. The polynomial coefficients have been adjusted to match these observations.
This is the code used in the JQuake software to calculate the position value:

public static double color2position(float hsv[]){    // Input : color in hsv space, float values between 0 and 1    double p = 0;    float h = hsv[0];    float s = hsv[1];    float v = hsv[2];    // Check if the color belongs to the scale    if (v > 0.1 && s > 0.75){        if (h > 0.1476){            p = 280.31*Math.pow(h,6) - 916.05*Math.pow(h,5) + 1142.6*Math.pow(h,4) - 709.95*Math.pow(h,3) + 234.65*Math.pow(h,2) - 40.27*h + 3.2217;        }        if (h <= 0.1476 && h > 0.001){            p = 151.4*Math.pow(h,4) -49.32*Math.pow(h,3) + 6.753*Math.pow(h,2) -2.481*h + 0.9033;        }        if (h <= 0.001){            p =   -0.005171*Math.pow(v,2) - 0.3282*v + 1.2236;        }    }    if (p < 0){        p = 0;    }    return p;}

To make people using other languages' lives easier (since other languages might return values from 0-255 instead of 0-1), I decided to write an equivalent C++ snippet. It takes R, G and B (0-255) values as arguments, converts it to 0-1 HSV values, and use the above algorithm to calculate the position value.

#include <math.h>#include <algorithm>#include <array>using namespace std;std::array<double, 3> getHsv(int r, int g, int b) {    double max = std::max(r, std::max(g, b));    double min = std::min(r, std::min(g, b));    if (min == max) {        return std::array<double, 3> {0, 0, max / 255};    }    double w = max - min;    double h = 0;    if (r == max) {        h = (g - b) / w;    }    if (g == max) {        h = ((b - r) / w) + 2;    }    if (b == max) {        h = ((r - g) / w) + 4;    }    if ((h *= 60) < 0) {        h += 360;    }    return std::array<double, 3> {h, (max - min) / max, max / 255};}double colorToPosition(int r, int g, int b) {    double p = 0;    std::array<double, 3> hsv = getHsv(r, g, b);    double h = hsv[0];    double s = hsv[1];    double v = hsv[2];    h /= 360;    if (v <= 0.1 || s <= 0.75) {        p = 0;    }    if (h > 0.1476) {        p = 280.31 * pow(h, 6) - 916.05 * pow(h, 5) + 1142.6 * pow(h, 4) - 709.95 * pow(h, 3) + 234.65 * pow(h, 2) - 40.27 * h + 3.2217;    } else if (h > 0.001) {        p = 151.4 * pow(h, 4) - 49.32 * pow(h, 3) + 6.753 * pow(h, 2) - 2.481 * h + 0.9033;    } else {        p = -0.005171 * pow(v, 2) - 0.3282 * v + 1.2236;    }    return std::max(0.0, p);}

Now that we have three interpolations, we can handle the position-to-number conversion.

Position to Number conversion

This step consists of converting the acquired position value p into a numerical value for various data types. The equations for doing so are as follows:

Seismic Intensity (Japan shindo scale decimal, Wikipedia)

Peak Ground Acceleration

Peak Ground Velocity

Peak Ground Displacement

Additional information on logarithmic scale

While creating conversion equations for logarithmic scales (PGA, PGV, ...) I was surprised to find that the ticks between 1, 2 and 5 are evenly spaced for every power of 10. I found that I couldn't figure out a mathematically valid way to project the values in such a consistent manner. I think the values are evenly spaced for display purposes only. Anyways, as you can see above, I used a decimal logarithmic scale.

If displayed correctly, the scale will look like this:

You can see that the space between 2 and 5 is larger than the space between 1 and 2. I know it isn't that important, but I'd like to respond to the thoughts of those who've tried to do so

Conclusion

In this article, I described how to derive numerical values from a GIF/PNG image from the seismic monitor website. To achieve this goal, I used polynomial interpolation. There are other methods that require less CPU, but are less robust against slight changes in color due to changes in the image format.
If you've read this far, thanks for reading this article, and have a nice day.

Additional resources

• NIED observation points database in JSON
• Official NIED website
• Download earthquake data from NIED