#! /usr/bin/env python
#
# Copyright (C) Apr 2011 Angelos Tzotsos <tzotsos@gmail.com>		  
# 
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.

try:
    from osgeo import gdal
    from osgeo.gdalconst import *
except ImportError:
    import gdal
    from gdalconst import *

try:
    import numpy
except ImportError:
    import Numeric as numpy

import sys

# =============================================================================
def Usage():
    print('Usage: gdal_accuracy.py classification_file truth_file')
    print('')
    sys.exit( 1 )

# =============================================================================
classfile = None
truthfile = None

# Parse command line arguments.
i = 1
while i < len(sys.argv):
    arg = sys.argv[i]

    if classfile is None:
        classfile = arg

    elif truthfile is None:
        truthfile = arg

    else:
        Usage()

    i = i + 1

if classfile is None:
    Usage()
if  truthfile is None:
    Usage()


classdataset = gdal.Open( classfile, GA_ReadOnly )
truthdataset = gdal.Open( truthfile, GA_ReadOnly )

if (classdataset.RasterCount != 1):
    print('Error: Classification file has more than 1 band')
    sys.exit( 1 )
if (truthdataset.RasterCount != 1):
    print('Error: Ground truth file has more than 1 band')
    sys.exit( 1 )
if ((classdataset.RasterXSize != truthdataset.RasterXSize) or
(classdataset.RasterYSize != truthdataset.RasterYSize)):
    print('Error: Classification and Ground truth file have different dimensions')
    sys.exit( 1 )
    

classband = classdataset.GetRasterBand(1)
truthband = truthdataset.GetRasterBand(1)
(classmin, classmax) = classband.ComputeRasterMinMax()
(truthmin, truthmax) = truthband.ComputeRasterMinMax()

if (truthmin != classmin) or (truthmax != classmax):
    print ('Error: Classification and Ground Truth values are incompatible')
    sys.exit(1)

TP = numpy.zeros((int(classmax)+1-int(classmin)))
FP = numpy.zeros((int(classmax)+1-int(classmin)))
FN = numpy.zeros((int(classmax)+1-int(classmin)))

for i in range(classband.YSize-1, -1, -1):
    scanline_cl = classband.ReadAsArray(0, i, classband.XSize, 1, classband.XSize, 1)
    scanline_tr = truthband.ReadAsArray(0, i, truthband.XSize, 1, truthband.XSize, 1)
    for c in range(0, classband.XSize):
        k=0
        for val in range(int(classmin), int(classmax)+1):
            if (scanline_cl[0][c] == val) and (scanline_tr[0][c] == val):
                TP[k] = TP[k] + 1
            elif (scanline_cl[0][c] == val) and (scanline_tr[0][c] != val):
                FP[k] = FP[k] + 1
            elif (scanline_cl[0][c] != val) and (scanline_tr[0][c] == val):
                FN[k] = FN[k] + 1
            k = k +1

Comp = TP/(TP+FN)
Correct = TP/(TP+FP)
Quality = TP/(TP+FP+FN)

print('-------------------------------------------------------------')
print('Classification Accuracy')
print('-------------------------------------------------------------')
print('  Class number  |  Completeness  |  Correctness  |  Quality  ')
for val in range(int(classmin), int(classmax)+1):
    print('       %d        |     %2.2f %%    |    %2.2f %%    |   %2.2f %%   ' %
          (val, 100*Comp[val], 100*Correct[val], 100*Quality[val]))
print('-------------------------------------------------------------')