如何使用 MATLAB 计算图像中曲线下的面积？

Question

alt text http://internationalpropertiesregistry.com/Server/showFile.php?file=%2FUpload%2Fstatistics.gifc49ca28823a561a41d09ef9adbb5e0c5.gif

x轴的单位是小时（h），一共是24小时。

y轴的单位是百万（m）。

如何以m*h为单位计算图片中红色曲线下的面积？

重要更新

只有图像可用（不是数据），我想以编程方式计算面积。

Answer 1

这是一个有趣的解决方案:)。顺便说一句，它使用需要一些用户交互的bwfill（类似于imfill）。

代码

%# Constants
gray_value_curve = 2;
gray_value_box = 3;
area_box_in_units = 10;

%# Read the image
I = imread('C:\p23\graph.gif');
%# Find the area of a unit block
figure(1);
imshow(I,[]);
[BS sq_elem] = bwfill;   
imshow(BS,[]);
%# Get the dimensions to make the estimate more accurate
X = zeros(size(BS));
X(sq_elem) = 1;
s = regionprops(X,'Area','BoundingBox');
block_area = s.Area + 2*(s.BoundingBox(3)-1) + 2*(s.BoundingBox(4)-1) + 4;

%#Find the area under the curve
I( ~(I == gray_value_curve | I == gray_value_box) ) = 0;
figure(2);
imshow(I,[]);
[BA area_curve_elem] = bwfill;
imshow(BA,[]);
%# Area under the curve
curve_area = numel(area_curve_elem);

%# Display the area in the required units
area = area_box_in_units*curve_area/block_area;
disp(area);

输出

113.5259

图 1 图2

Answer 2

创建全自动解决方案的困难在于，它需要您将有关要处理的输入图像的某些假设硬编码到您的解决方案中。如果这些假设不适用于所有您可能遇到的潜在图像，那么全自动解决方案将不会给出值得信赖的结果，并尝试扩展全自动解决方案以处理 所有可能的输入都可能导致它膨胀成难以理解的复杂代码。

如果您对输入图像特征的可变性存有疑问，具有一些用户交互的Jacob's 之类的解决方案通常是最好的。如果您可以确定输入图像的特征遵循一组严格的规则，那么可以考虑使用自动化解决方案。

例如，下面是我编写的一些自动化代码，用于近似图表中红色曲线下方的区域。由于我使用上图作为指导，因此它必须满足许多条件才能工作：

• 绘制线的红色像素必须在图像中唯一地描述为包含等于 0 的绿色和蓝色分量以及等于 1 的红色分量。
• 网格线的绿色像素必须在图像中被唯一地描述为包含小于 1 的红色和蓝色分量以及等于 1 的绿色分量。
• 轴线的蓝色像素必须在图像中被唯一地描述为包含等于 0 的红色和绿色分量以及等于 1 的蓝色分量。
• 网格线和轴线必须始终在水平或垂直方向上完全对齐。
• 网格线的长度必须跨越图片宽度和高度的一半以上。
• x 轴必须是图像中最长的水平蓝线。
• 网格线的宽度必须始终为 1 像素。

在输入图像上满足以上条件，无需用户输入，可以使用下面的代码来逼近红色曲线下的区域：

[img,map] = imread('original_chart.gif');  %# Read the indexed image
[r,c] = size(img);                         %# Get the image size

redIndex = find((map(:,1) == 1) & ...    %# Find the red index value
                (map(:,2) == 0) & ...
                (map(:,3) == 0))-1;
greenIndex = find((map(:,1) < 1) & ...   %# Find the green index value
                  (map(:,2) == 1) & ...
                  (map(:,3) < 1))-1;
blueIndex = find((map(:,1) == 0) & ...   %# Find the blue index value
                 (map(:,2) == 0) & ...
                 (map(:,3) == 1))-1;

redLine = (img == redIndex);      %# A binary image to locate the red line
greenLine = (img == greenIndex);  %# A binary image to locate the grid lines
blueLine = (img == blueIndex);    %# A binary image to locate the axes lines

w = mean(diff(find(sum(greenLine,1) > r/2)));  %# Compute unit square width
h = mean(diff(find(sum(greenLine,2) > c/2)));  %# Compute unit square height
squareArea = w*h;                              %# Compute unit square area

[maxValue,maxIndex] = max(redLine);          %# Find top edge of red line
x = find(maxValue > 0);                      %# Find x coordinates of red line
y = maxIndex(maxValue > 0);                  %# Find y coordinates of red line
[maxValue,maxIndex] = max(sum(blueLine,2));  %# Find row index of x axis
y = maxIndex-y;                              %# Zero the y coordinate
totalArea = trapz(x,y)/squareArea;           %# Compute the area under the curve

结果如下：

squareArea = 460.6 square pixels
totalArea = 169.35 m*h

解释：

我将详细说明计算w所涉及的步骤：

1. 二进制图像greenLine 使用函数SUM 沿每一列求和，得到一个1-by-c 向量，其中每个元素是图像每列中有多少网格线像素的计数。
2. 此向量中大于r/2（图像中行数的一半）的元素表示图像中包含垂直网格线的列。这些列的索引是使用函数FIND 找到的。
3. 使用函数DIFF 可以找到这些列索引之间的成对差异。这给出了一个向量，其中包含网格线之间空间的宽度（以像素为单位）。
4. 最后，函数MEAN 用于计算图像中所有网格线之间空间的平均宽度。

在计算h 时，唯一的区别是对每一行进行求和，并将r/2 替换为c/2（图像中列数的一半）。

Answer 3

由于您只有可用的图像，我建议您通过肉眼进行积分：计算红线下方的整个正方形的数量。

对于红线相交的每个方格，根据红线下方的数量决定是否将其包括在计数中。不要试图估计有多少正方形位于红线下方，这充其量只会让您产生更准确的错觉。

编辑：我为你数了绿色方块，答案是 168 m.h

Answer 4

因为这看起来不像是一个可以积分的“函数”，所以我会使用数值积分技术。我总是偏爱trapz，它使用“trapezoidal rule”进行数值积分。

类似：

area = trapz(data);

应该足够了。

希望对你有帮助，

会