mAP以及AP计算方式

1、前言

　　在目标检测中我们常常使用AP（Average Precision）作为模型对某种目标精测精度的评价指标，该如何计算AP呢？

　　AP就是P—R曲线下的面积，我们需要做的就是根据不同的置信度阈值(p_threshold)，计算出这模型得到的预测框的（R，P），然后作出P—R曲线，并求解面积，就能得到目标检测模型对该检测种类的AP。(在VOC2010之后，计算AP需要对做出来的P—R曲线做一个平滑，之后会提到)

　　和分类模型计算P（精确度）R（召回度）不一样的是（分类模型计算P，R可以参看我这篇博文AUC、精确率和召回率），在目标检测中我们是没法知道我们模型出来的预测框的真实值是什么的(甚至标记的真实框的个数和预测框的个数都是不匹配的)，该怎样计算P和R呢？其实计算P，R最重要的指标只有3个，分别是：

TP	真实值为正——>预测值为正
FP	真实值为负——>预测值为正
FN	真实值为正——>预测值为负

　　因此在目标检测中我们只需要计算出以上三个指标就可以了。

2、例子

　　我们直接用例子来进行解释，假设我们有5个GT（Ground_Truth就是目标检测中人工标注的对象的矩形框，可以认为是真实物体的标签，下面用BB1—BB5来表示），我们模型对该种识别种类返回了10个预测框，接下来我们计算该种类的AP：

　　1、根据预测框计算与5个GT的IOU，如果IOU大于iou_threshold(IOU阈值，不要和前面的p_threshold弄混了)则标记为1，否则标记为0。比如我们预测出来10个框，根据iou_threshold可以得到这些框的IOU结果分别为（1，0，1，0，1，1，1，0，0，1），还有这些框的置信度分别为[0.9，0.85，0.7，0.6，0.45，0.25，0.2，0.15，0.13，0.12]（需要对框从小到大排序），并且我们还知道这些预测框和哪一个GT在做IOU，[BB1，BB2，BB1，BB2，BB3，BB4，BB3，BB1，BB4，BB5]，如果做一个表来表示结果应该如下图：

编号	GT	置信度	IOU结果
1	BB1	0.9	1
2	BB2	0.85	0
3	BB1	0.7	1
4	BB2	0.6	0
5	BB3	0.45	1
6	BB4	0.25	1
7	BB3	0.2	1
8	BB1	0.15	0
9	BB4	0.13	0
10	BB5	0.12	1

　　2、要注意虽然预测出来有10个框，但是我们实际上只有5个GT（也就是只有5个物体），换句话说其实你预测的10个框里可能有对同一个GT的重复预测（可能有几个预测框都和同一个GT满足IOU>iou_threshold，比如1号框和3号框都是同时对BB1的预测）。

　　3、接下来我们就是需要根据不同的p_threshold(置信度阈值)，计算出这10个预测框的（R，P）点，然后作出P——R曲线，并求解面积。问题就在于如何求对应p_threshold下的P和R，下面我以p_threshold=0.6为例，计算上面表格中10个框的（R，P）具体做法。

由于编号1，2，3，4的预测框的置信度都大于等于p_threshold，因此我同意他们的IOU结果(即认为编号1预测框的标签是BB1的准确框，编号2不是BB2的准确框，编号3也是BB1的准确框，编号4不是BB2的准确框)，并且我认为前面四个框(编号1-4)都预测为正例。
其余的预测框我都不认同他们的IOU结果(即我不认为编号5的预测框是BB3的准确框，我也不认为编号8不是BB1的准确框)，并且我仍然把剩下的6个框都预测为负例。
由上面两种看法，我们就可以得到当p_threshold=0.6时：这样我们就可以算出P=TP/(TP+FP)，R=TP/(TP+FN)
- TP = 1（只有编号1是正确预测，编号3虽然也是正确预测但是属于重复预测了归纳到FP中，谁让它置信度小呢）
- FP = 3（因为编号2的IOU结果表示它为负类，但是我预测为正类，因此编号2是FP，同理编号4也是。编号3是认为重复预测了，因为编号1已经正确预测了BB1）
- FN = 4（这个指标可以直接由GT个数-TP个数得到，因为你可以理解为真实的有5个正例，但是你只预测对了1个(TP=1)，所以剩下的框都预测错了，因此FN=GT-TP）

　　上面的例子说完其实可能有疑问说为什么不取分析编号5-10的框呢？置信度小于我设定的阈值p_threshold，我没有办法认同他们的IOU结果因此我不知道他们的IOU结果是对的还是错的，因此我们并不能也没有必要去分析他们。（因为我们并不在意TN是多少）

　　另一个疑问是，为什么你要把置信度大于p_threshold的框都预测为正例，其余的是负例，难道不是应该尊重IOU结果么？我的想法是模型既然输出的这些预测框，那么肯定是模型认为这些预测框都是目标框，但是置信度小于了我设定的p_threshold自然是应该认为模型识别为负例了。

　　其实，以上的看法都是我我自己根据网上博文得到的结果的理解，并且计算之后发现这一套逻辑是符合最后结果的，为了方便大家及以和自己去理解这个指标。

3、AP的计算

　　通过上面的例子，我们知道不同的p_threshold可以计算得出不同的（R，P），由这些点可以画出图像(这个图像就是P—R图像)，我们可以知道P会随着R上升而下降，根据VOC2010的计算公式，我们需要对P—R曲线做一个平滑，具体的方式很简单。就是每一个点的R1所对应的P1调整为 $P1 = \max\limits_{R>R1}P(R)$ 。用图像来表示就是从图一平滑成图二。

图一

图二

　　根据平滑后的P—R图像的结果我们就可以很轻松的计算面积了(其实就是矩形求面积)。

4、mAP的计算

　　对于mAP（mean of Average Precision）的计算，其实就是在目标检测中我们可能是在检测多个目标，因此我们对每种类型的目标都可以计算出这个种类的AP，最后再对各个种类进行求平均就可以了。

5、代码及解释

　　根据上面的分析我们可以知道，对我们在对P-R曲线做光滑的时候其实就是在保证P随着R的增加应该是单调不增的函数，如何得到呢？我们可以对precision列表(根据recall从小到大对应排序)采用从末尾到首端不断取最大值就可以了。（也就是完整代码代码的第48行—第53行）

　　其实我们还可以发现，我们并不需要对所有的p_threshold从0到1的所有取值去得到（R,P）点，因为只有当p_threshold大于某个框的置信度的时候，P和R才会改变，因此我们算(R,P)点只需要计算p_threshold等于各个框的置信度的时候就可以了。（并且我们会发现随着p_threshold的增加，racall是一定增加的。）

　　在完整代码的第48行到第57行我们求面积的时候运用了错位的方式巧妙的找到了precision突变的那些位置，为了更好的展示这部分的工作，我虚构了一些precision和recall来让大家更好的理解（可以结合图二来理解，数据不同）。代码如下：

mpre = np.array([0.,0.98,0.88,0.88,0.25,0.3,0.32,0.21,0.33,0.1,0.])
for i in range(mpre.size - 1, 0, -1):
    mpre[i - 1] = np.maximum(mpre[i - 1], mpre[i])
print(mpre)

[0.98 0.98 0.88 0.88 0.33 0.33 0.33 0.33 0.33 0.1  0.  ]

mrec = np.array([0.,0.12,0.12,0.56,0.56,0.56,0.75,0.75,0.82,0.93,1.0])
i = np.where(mrec[1:] != mrec[:-1])[0]
print(i)
print('mrec[i + 1]:',mrec[i + 1])
print("mrec[i]:",mrec[i])
print("mpre[i + 1]:",mpre[i + 1])

[0 2 5 7 8 9]
mrec[i + 1]: [0.12 0.56 0.75 0.82 0.93 1.  ]
mrec[i]: [0.   0.12 0.56 0.75 0.82 0.93]
mpre[i + 1]: [0.98 0.88 0.33 0.33 0.1  0.  ]

　　上面的示例可以看出来完整代码的第57行到第60行是可以完全符合我们求取矩形面积框的需求的。

　　因为代码是基于VOC数据集编写的，所以我们最好了解一下VOC数据集的XML文件格式：

<annotation>
    <folder>VOC2012</folder>                           
    <filename>2007_000392.jpg</filename>              //文件名
    <source>                                         //图像来源（不重要）
        <database>The VOC2007 Database</database>
        <annotation>PASCAL VOC2007</annotation>
        <image>flickr</image>
    </source>
    <size>                          //图像尺寸（长宽以及通道数）                        
        <width>500</width>
        <height>332</height>
        <depth>3</depth>
    </size>
    <segmented>1</segmented>        //是否用于分割（在图像物体识别中01无所谓）
    <object>                        //检测到的物体
        <name>horse</name>          //物体类别
        <pose>Right</pose>          //拍摄角度
        <truncated>0</truncated>    //是否被截断（0表示完整）
        <difficult>0</difficult>    //目标是否难以识别（0表示容易识别）
        <bndbox>                    //bounding-box（包含左下角和右上角xy坐标）
            <xmin>100</xmin>
            <ymin>96</ymin>
            <xmax>355</xmax>
            <ymax>324</ymax>
        </bndbox>
    </object>
    <object>                        //检测到多个物体
        <name>person</name>
        <pose>Unspecified</pose>
        <truncated>0</truncated>
        <difficult>0</difficult>
        <bndbox>
            <xmin>198</xmin>
            <ymin>58</ymin>
            <xmax>286</xmax>
            <ymax>197</ymax>
        </bndbox>
    </object>
</annotation>

　　完整代码如下：

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

import numpy as np
import os
import pickle
import xml.etree.ElementTree as ET
 
def parse_rec(filename):
    """
    解析PASCAL VOC xml文件
    return：[{'name': xxx, 'bbox': [xmin, ymin, xmax, ymax]},{},....]
    """
    tree = ET.parse(filename)
    objects = []
    for obj in tree.findall('object'):
        obj_struct = {}
        obj_struct['name'] = obj.find('name').text
        obj_struct['pose'] = obj.find('pose').text
        obj_struct['truncated'] = int(obj.find('truncated').text)
        obj_struct['difficult'] = int(obj.find('difficult').text)
        bbox = obj.find('bndbox')
        obj_struct['bbox'] = [int(bbox.find('xmin').text),
                              int(bbox.find('ymin').text),
                              int(bbox.find('xmax').text),
                              int(bbox.find('ymax').text)]
        objects.append(obj_struct)
  
    return objects
 
 
def voc_ap(rec, prec, use_07_metric=False):
    """
    给定recall和precision之后计算返回AP，其中recall是从小到大排序，precision每一个元素是对应排序后的recall的值
    rec: np.array
    prec: np.array
    use_07_metric:如果为True则采用07年的方式计算AP
    """
    if use_07_metric:  # VOC在2010之后换了评价方法，所以决定是否用07年的
        ap = 0.
        for t in np.arange(0., 1.1, 0.1):  #  07年的采用11个点平分recall来计算
            if np.sum(rec >= t) == 0:
                p = 0
            else:
                p = np.max(prec[rec >= t])  # 取一个recall阈值t之后最大的precision
            ap = ap + p / 11.  # 将11个precision加和平均
    else:  
        # 这里是使用VOC2010年后的方法求mAP，计算光滑后PR曲线的面积，不再是固定的11个点
        # 在rec的首和尾添加值来完成更好完成
        # 在prec的尾部添加0为了更好得到“光滑”后的prec的值
        mrec = np.concatenate(([0.], rec, [1.]))  # recall和precision前后分别加了一个值，因为recall最后是1，所以
        mpre = np.concatenate(([0.], prec, [0.])) # 右边加了1，precision加的是0
 
        # 调整mpre，从后往前取最大值，保证prec单调不增。
        for i in range(mpre.size - 1, 0, -1):
            mpre[i - 1] = np.maximum(mpre[i - 1], mpre[i])  
 
        # 计算PR曲线下的面积
        # X轴为R（recall的值）
        i = np.where(mrec[1:] != mrec[:-1])[0]  # 返回了所有改变了recall的点的位置
        # 求每个矩形的面积和
        # 具体理解见前文解释
        ap = np.sum((mrec[i + 1] - mrec[i]) * mpre[i + 1])
    return ap
 
# 计算每个类别对应的AP，mAP是所有类别AP的平均值
# 主要是处理得到rec, prec数组
def voc_eval(detpath,
             annopath,
             imagesetfile,
             classname,
             cachedir,
             ovthresh=0.5,
             use_07_metric=False):
    """
    return：rec, prec, ap
    Top level function that does the PASCAL VOC evaluation.
    detpath: 检测结果的文件路径。检测结果文件的每一行应该是：img_ID，置信度，xmin, ymin, xmax, ymax
            detpath应该是这样的字符串'./results/comp4_det_test_{}.txt'
            detpath.format(classname) should produce the detection results file.
    annopath: Path to annotations
            annopath应该是这样的字符串"dataset/voc/VOC2007/Annotations/{}.xml"
            annopath.format(imagename) should be the xml annotations file.
    imagesetfile: 储存了图片名的text，每一行是一张图片的名。'dataset/voc/VOC2007/ImageSets/Main/test.txt'
    classname: 类别名
    cachedir: 用于存储注解(annotations)的路径，生成一个pickle_file
    [ovthresh]: IOU_threshold (default = 0.5)
    [use_07_metric]: 是否使用07年的计算ap的方式(default False)
          
    """
    # 第一步获得各图片的GT
    # 如果不存在注释路径的文件夹，则先创建文件夹
    if not os.path.isdir(cachedir):
        os.mkdir(cachedir)
    cachefile = os.path.join(cachedir, 'annoys.pkl')
  
    # 读取图片名，并储存在列表中
    with open(imagesetfile, 'r') as f:
        lines = f.readlines()
    imagenames = [x.strip() for x in lines]
  
    if not os.path.isfile(cachefile):
        # 如果没有现成的解析好的图片GT数据则需要从xml文件中解析出来，并保存成pkl文件
        # 这里提取的是所有测试图片中的所有object gt信息, 07年的test真实标注是可获得的，12年就没有了
        recs = {}
        for i, imagename in enumerate(imagenames):
            recs[imagename] = parse_rec(annopath.format(imagename)) # 获取图片中对应的GT的解析
            if i % 100 == 0:
                print('Reading annotation for {:d}/{:d}'.format(
                    i + 1, len(imagenames)))
        # save
        print('Saving cached annotations to {:s}'.format(cachefile))
        with open(cachefile, 'wb') as f:
            pickle.dump(recs, f)
    else:
        # 如果有现成图片GT数据，则直接读取
        with open(cachefile, 'rb') as f:
            recs = pickle.load(f)
  
    # 从上面的recs中提取出我们要判断的那类目标的标注信息(GT)
    class_recs = {}
    npos = 0
    for imagename in imagenames:
        R = [obj for obj in recs[imagename] if obj['name'] == classname]    # [obj，obj，....] 每个obj={'name': xxx, 'bbox': [xmin, ymin, xmax, ymax]}
        bbox = np.array([x['bbox'] for x in R])     # 二维数组，(number_obj,4),该张图片有number_obj个类别为classname的目标框
        difficult = np.array([x['difficult'] for x in R]).astype(np.bool)
        det = [False] * len(R)        # 该图片中该类别对应的所有bbox的是否已被匹配的标志位
        npos = npos + sum(~difficult) # 累计所有图片中的该类别目标的GT总数，不算diffcult
  
        class_recs[imagename] = {'bbox': bbox,
                                 'difficult': difficult,
                                 'det': det}
  
    # 第二步读取模型识别的结果
    detfile = detpath.format(classname)
    # 读取相应类别的检测结果文件，每一行对应一个检测目标
    with open(detfile, 'r') as f:
        lines = f.readlines()   # 读取所有行
  
    splitlines = [x.strip().split(' ') for x in lines]  # 处理为[[image_id, 置信度, xmin, ymin, xmax, ymax ],[]...]
    image_ids = [x[0] for x in splitlines]
    confidence = np.array([float(x[1]) for x in splitlines])        # 一维数组
    BB = np.array([[float(z) for z in x[2:]] for x in splitlines])  # 二维数组，（number_bbox,4）
  
    # sort by confidence 按置信度由大到小排序
    sorted_ind = np.argsort(-confidence)    # 获得Indx
    # sorted_scores = np.sort(-confidence)
    BB = BB[sorted_ind, :]  # 对BB重排序
    image_ids = [image_ids[x] for x in sorted_ind]  # 对image_ids重排序
  
    # 记下dets并对每个image打上标注是TP还是FP
    nd = len(image_ids) # 检测结果文件的行数
    tp = np.zeros(nd)   # 用于标记每个检测结果是tp还是fp
    fp = np.zeros(nd)
    for d in range(nd):
       # 取出该条检测结果所属图片中的所有ground truth
        R = class_recs[image_ids[d]]    # 其实image_id就是image_name，R={'bbox': bbox(二维数组),'difficult': difficult,'det': [bool]}
        bb = BB[d, :].astype(float)     # bb一维数组
        ovmax = -np.inf
        BBGT = R['bbox'].astype(float)  # 二维数组
  
        if BBGT.size > 0:
            # compute overlaps  计算与该图片中所有ground truth的最大重叠度(IOU)
            # intersection
            ixmin = np.maximum(BBGT[:, 0], bb[0])   # 一维
            iymin = np.maximum(BBGT[:, 1], bb[1])
            ixmax = np.minimum(BBGT[:, 2], bb[2])
            iymax = np.minimum(BBGT[:, 3], bb[3])
            iw = np.maximum(ixmax - ixmin + 1., 0.)
            ih = np.maximum(iymax - iymin + 1., 0.)
            inters = iw * ih    # 一维
  
            # 重叠部分面积一维
            uni = ((bb[2] - bb[0] + 1.) * (bb[3] - bb[1] + 1.) +
                   (BBGT[:, 2] - BBGT[:, 0] + 1.) *
                   (BBGT[:, 3] - BBGT[:, 1] + 1.) - inters)
  
            overlaps = inters / uni     # 计算得到检测结果的这个框与该张图片的所有该类比的GT的IOU，一维
            ovmax = np.max(overlaps)
            jmax = np.argmax(overlaps)
  
        # 这里就是具体的分配TP和FP的规则了
        if ovmax > ovthresh:  # 如果最大的重叠度大于一定的阈值（IOU_threshold）
            if not R['difficult'][jmax]: # 如果最大重叠度对应的ground truth为difficult就忽略，因为上面npos就没算
                                                              
                if not R['det'][jmax]: # 如果对应的最大重叠度的ground truth以前没被匹配过则匹配成功，即tp
                    tp[d] = 1.
                    R['det'][jmax] = 1  # 表示框被匹配过了
                else:  # 若之前有置信度更高的检测结果匹配过这个ground truth，则此次检测结果为fp
                    fp[d] = 1.
        else:
            # 该图片中没有对应类别的目标ground truth或者与所有ground truth重叠度都小于阈值
            fp[d] = 1.
  
    # 计算 precision recall
    fp = np.cumsum(fp) # 累加函数np.cumsum([1, 2, 3, 4]) -> [1, 3, 6, 10]
    tp = np.cumsum(tp)
    rec = tp / float(npos)  #  tp/GT,也就得到了voc_ap函数所需要的rec了
    # avoid divide by zero in case the first detection matches a difficult
    # 避免除以零
    prec = tp / np.maximum(tp + fp, np.finfo(np.float64).eps)
    ap = voc_ap(rec, prec, use_07_metric)
  
    return rec, prec, ap

　　下面给出函数voc_eval的具体输入格式以及需要的各个文件的格式：

detpath = '/home/g4/桌面/detection_result_{}.txt'
annopath = '/home/g4/桌面/labels/{}.xml'
imagesetfile = '/home/g4/imagesetfile.txt'
classname = '烟雾'
cachedir = '.'
rec, prec, ap = voc_eval(detpath,annopath,imagesetfile,classname,cachedir)
ap