Hi,
Markus Neteler wrote:
On Fri, Sep 11, 2009 at 5:25 PM, Georg Kaspar <[email protected]> wrote:
Dear GRASS developers,
Today I finished a first attempt at writing a GRASS script in python for
probabilistic label relaxation.
Perhaps you could post an example here (ideally based on the North Carolina
data set) to facilitate testing?
here:
http://www.file-upload.net/download-1894443/plr_in_grass.pdf.html
you will find a small documentation including examples from the NC data set.
I also attached an updated version of the code which contained some
nasty bugs - apart from being slow it should work now ;)
best regards,
Georg
#!/usr/bin/env python
#
#############################################################################
#
# MODULE: i.plr.py
# AUTHOR(S): Georg Kaspar
# PURPOSE: Probabilistic label relaxation, postclassification filter
# COPYRIGHT: (C) 2009
#
# This program is free software under the GNU General Public
# License (>=v2). Read the file COPYING that comes with GRASS
# for details.
#
#############################################################################
#%Module
#% description: Probabilistic label relaxation, postclassification filter
#%End
#%option
#% key: group
#% type: string
#% description: image group to be used
#% required : yes
#%end
#%option
#% key: subgroup
#% type: string
#% description: image subgroup to be used
#% required : yes
#%end
#%option
#% key: sigfile
#% type: string
#% description: Path to sigfile
#% required : yes
#%end
#%option
#% key: output
#% type: string
#% description: Name for resulting raster file
#% required : yes
#%end
#%option
#% key: iterations
#% type: integer
#% description: Number of iterations (1 by default)
#% required : no
#%end
#%option
#% key: ntype
#% type: integer
#% description: type of neighbourhood (4(default) or 8)
#% required : no
#%end
import sys
import os
import numpy
import grass.script as grass
from osgeo import gdal, gdalnumeric, gdal_array
from osgeo.gdalconst import GDT_Byte
def getMetadata():
env = grass.gisenv()
global GISDBASE
global LOCATION_NAME
global MAPSET
GISDBASE = env['GISDBASE']
LOCATION_NAME = env['LOCATION_NAME']
MAPSET = env['MAPSET']
def splitSignatures(path, sigfile):
# split signature file into subfiles with 1 signature each
stream_in = open(path + sigfile, "r")
stream_in.next() # skip first line
counter = 0
stream_out = open(path + "plr_foo.sig", "w")
for line in stream_in:
if line[0] == "#":
stream_out.close()
counter += 1
stream_out = open(path + "plr_" + str(counter) + ".sig", "w")
stream_out.write("#produced by i.plr\n")
stream_out.write(line)
stream_out.close()
stream_in.close()
return counter
def normalizeProbabilities(counter):
arg = ""
for i in range(counter):
arg = arg + "+plr_rej_" + str(i)
arg = arg.strip('+')
print "calculating multiplicands, arg=" + arg
grass.run_command("r.mapcalc", multiplicand = "1./(" + arg + ")")
for i in range(counter):
print "normalizing probabilities for class " + str(i)
grass.run_command("r.mapcalc", plr_rej_norm = "plr_rej_" + str(i) + "*multiplicand")
grass.run_command("g.rename", rast="plr_rej_norm,plr_rej_norm_" + str(i))
def getProbability(a,b):
# TODO: Implement this!!!
if a == b:
return 1
else:
return 0.5
def cleanUp(path):
os.system("rm " + path + "/plr_*.*")
os.system("rm /tmp/plr_*.*")
grass.run_command("g.mremove", flags="f", quiet=True, rast="plr_*")
def plr_filter(probabilities, width, height, classes, type):
# create an empty n-dimesional array containing results
results = numpy.ones((classes,height,width))
print "starting label relaxation"
progress = 0
# for each pixel (except border)
for y in range(height):
p = int(float(y)/height * 100)
if p > progress:
print '\r' + str(p) + '%'
progress = p
for x in range(width):
# for all classes create neighbourhood and extract probabilities
for c in range(1, classes+1):
if (x == 0) or (x == width-1) or (y == 0) or (y == height-1):
results[c-1,y,x] = probabilities[c-1,y,x]
else:
if type == 8:
q = neighbourhoodFunction8(probabilities, x, y, c, classes)
else:
q = neighbourhoodFunction4(probabilities, x, y, c, classes)
p = probabilities[c-1, y, x]
# resulting cell contains the product of class probability and
# neighbourhood-function
results[c-1,y,x] = p * q
print ""
return results
def neighbourhoodFunction4(probabilities, x, y, z, classes):
n = []
neighbours = [[x-1,y],[x,y],[x+1,y],[x,y-1],[x,y+1]]
for i in neighbours:
l = []
# for each possible class
for c in range(1, classes+1):
l.append(getProbability(z, c) * float(probabilities[c-1,i[1],i[0]]))
n.append(sum(l))
return sum(n)
def neighbourhoodFunction8(probabilities, x, y, z, classes):
n = []
for j in range(y-1, y+2):
for i in range(x-1, x+2):
l = []
# for each possible class
for c in range(1, classes+1):
l.append(getProbability(z, c) * float(probabilities[c-1,j,i]))
n.append(sum(l))
return sum(n)
def createMap(probabilities, width, height, classes):
print "retrieving class labels"
results = numpy.ones((height,width))
progress = 0
for y in range(height):
p = int(float(y)/height * 100)
if p > progress:
print '\r' + str(p) + '%'
progress = p
for x in range(width):
max_class = 1
max_val = probabilities[0,y,x]
for c in range(2, classes+1):
current_val = probabilities[c-1,y,x]
if current_val > max_val:
max_val = current_val
max_class = c
results[y,x] = max_class
#results = probabilities.max(0)
return results;
def export(array, trans, proj):
driver = gdal.GetDriverByName('GTiff')
out = driver.Create('/tmp/plr_results.tif', array.shape[1], array.shape[0], 1, GDT_Byte)
out.SetGeoTransform(trans)
out.SetProjection(proj)
gdal_array.BandWriteArray(out.GetRasterBand(1), array)
def main():
# fetch parameters
group = options['group']
subgroup = options['subgroup']
sigfile = options['sigfile']
output = options['output']
iterations = options['iterations']
ntype = options['ntype']
if iterations == "":
iterations = 1
iterations = int(iterations)
if ntype == "":
ntype = 4
ntype = int(ntype)
# fetch Metadata
getMetadata()
# split sigfiles
sigpath = GISDBASE + "/" + LOCATION_NAME + "/" + MAPSET + "/group/" + group + "/subgroup/" + subgroup + "/sig/"
counter = splitSignatures(sigpath, sigfile)
print "found " + str(counter) + " signatures"
l = []
for i in range(1, counter+1):
# extract probabilities
print "extracting probabilities for class " + str(i)
grass.run_command("i.maxlik", group=group, subgroup=subgroup, sigfile="plr_" + str(i) + ".sig", clas="plr_class" + str(i), reject="plr_rej_" + str(i))
# export from GRASS
print "exporting probabilities for class " + str(i) + " to /tmp"
grass.run_command("r.out.gdal", inp="plr_rej_" + str(i), out="/tmp/plr_rej_" + str(i) + ".tif")
# import via gdal
print "reading file"
tif = gdal.Open("/tmp/plr_rej_" + str(i) + ".tif")
l.append(tif.ReadAsArray())
if i == 1:
width = l[0].shape[1]
height = l[0].shape[0]
trans = tif.GetGeoTransform()
proj = tif.GetProjection()
# create n-dimensional array
print "creating 3D-array"
probabilities = numpy.array(l)
print "Image size: " + str(width) + "x" + str(height)
print "using " + str(ntype) + "-neighbourhood"
# invoke relaxation process
results = probabilities.copy()
for i in range(int(iterations)):
print str(int(iterations)-i) + " iteration(s) to go..."
results = plr_filter(results, width, height, counter, ntype)
labels = createMap(results, width, height, counter)
# exporting results
print "exporting results to /tmp"
export(labels, trans, proj)
# import via gdal into GRASS
print "reading results"
grass.run_command("r.in.gdal", inp="/tmp/plr_results.tif", out=output)
# clean up
print "removing temporary files"
cleanUp(sigpath)
if __name__ == "__main__":
options, flags = grass.parser()
main()
_______________________________________________
grass-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev
