'''
Mean Filter

Created by:  Nathan Jennings
Updated on:  November 4, 2007

Mean Filter using Numeric arrays, Gdal built for Windows, Python 2.4

This filter uses a kernel structure to process pixel elements stored in
a Numeric array.

Input - TIF file
Output - TIF file as described in Gdal
kernel size - 3x3 (3), 5x5 (5), etc

NOTES:  This function will only work on Windows OS and uses the Python 2.4
gdal build for Python 2.4 as found on www.gdal.org, Frank Warmerdam. Numeric
Python (Numeric) is an older version of Numpy (or SciPy) that can be used
with this version of gdal 1.4.2.

I developed this function using a 6 x 400 x 600 Landsat TM subset.
It took ~ 5 minutes, which I think is long.  I am not using a look up table
or any other optimization, but I have read that Python processes slower than
C++ or other lower level programming language.

This same structure can be used to put in any other kind of kernel based processing
such as median, guassian, or custom kernel.  Most of the code will likely go inside
the inner most "for loops."  Other intermediate calculations may be needed and will
need to be placed accordingly.

I will look into optimization routines to incorporate or if others have suggestions
you can contact me at nate@deltateck.com.  I would like to hear about them or being
able to do something similar in C++, provided gdal will work with it.  I would be
using the free C++ compiler (g++ or similar) and not writing the actual bindings
to build gdal.

'''

import gdal
import Numeric
import os, sys
import time


#Open Multi-band image"


print "Starting Mean Filter..."

img = gdal.Open(r'g:\\remote_sensing_images\\Landsat\\TM\\sactmclp_sub.tif')


#Load Image into a Numeric (Python) Array
#The line below is necessary because I use it below.
#You don't have to use an array for this, but you can use various
#gdal objects to get the bands, rows, and columns of an image.  See below.
#dal requires processing by band....see below.

a = img.ReadAsArray()


#Assign an output file.  I have it hard coded here.

outfile = "g:\\remote_sensing_images\\Landsat\\TM\\test_mean_img.tif"

#Create output file parameters

#The value here is one of many values for different formats of imagery
#check out the website above for more details.

format = "GTiff"
out_driver = gdal.GetDriverByName( format )

#Create Copy is a gdal function that basically makes a copy of the input image
#This will also keep any georeferencing with the image.
#If Create is used, then you will loose the georeference and have to add
#it in either programmatically or through other software.

outputimg = out_driver.CreateCopy(outfile, img, 0)


#The following shows how to report out the bands, rows columns
#Once it is in an array.
#You can also report out below using gdal objects such as
#RasterCount, RasterX, and RasterY (i.e. img.RasterCount)

#The following is not necessary to run the function, but for illustrative
#purposes and to show some of my own investigation.

print a.shape   #Prints the dimensions of the iamge array
print "Num Bands: " + str(a.shape[0])
print "Num Rows: " + str(a.shape[1])
print "Num Cols: " + str(a.shape[2])

#Like above, the following variables could have been assigned from using
#the RasterCount, RasterX, and RasterY from the input image

nbnds = a.shape[0]
nrows = a.shape[1]
ncols = a.shape[2]


#Filter (kernel) size - 3x3 - filter = 3
filter = 3
Window_Total = filter*filter
offset = (filter-1)/2

#This sets up the proper row and column to start the processing
#Since the method operates by using a kernel, the index for the
#rows and columns need to be offset, which is determined by the kernel size

row = offset   #For a 3x3 kernel, the offset will be 1; 5x5 will be 2, etc.
col = offset

bnd = 1        #This initializes the band for loop


#Just a counter.  I wanted to see how many calculations were being done, approximately
CalcCount = 0

#The following loops over the bands, row,s and columns, calculates the mean
#and assigns them to the respective output

#Need to include img.RasterCount + 1 because the range option for the for loop
#loops up to, but not including the end range (img.RasterCount)

for bnd in range(bnd, img.RasterCount+1):

    #Gdal operates on a given band...this gets the band
    #Then stores it in a Numeric array, so that it can be processed

    in_band = img.GetRasterBand(bnd)
    band_array = in_band.ReadAsArray()

    #gets the band of the output
    #and allocates an array for it.  The array is where the results of the
    #mean will be stored
    
    out_band = outputimg.GetRasterBand(bnd)  
    out_array = out_band.ReadAsArray()
    


    row = offset            #Resets the row for the next band

    #After reading in a band to an array, loop over rows and columns

    for row in range(row, nrows-offset):

        col = offset        #Resets the column to the first column of the next row
        
        for col in range(col, ncols-offset):

            #Work on the kernel of pixels within a given band
            #This looping gets to a target pixel to do calculations

            winrow = row - offset   #Set the initial row for the kernel
            
            WinSum = 0
            
            CalcMean = 0    #Initialize the CaclMean variable
            
            for winrow in range(winrow, row+offset+1):
                wincol = col - offset    #Make sure to start the next row at column 0
                for wincol in range(wincol, col+offset+1):
                    

                    #Sum the pixels in the kernel
                    #The array stores a specific value in the array (i.e. the pixel value)
                    WinSum += band_array[winrow][wincol]
                    
                    
            #Notice the CalcMean is indented to calculate within the "col for loop"
            #not within the "wincol loop"
            
            CalcMean = WinSum / Window_Total
            

            #Assign mean within kernel to output array
            #This is offset by 1 row and 1 column as a result of the kernel size being used
            
            out_array[row][col] = CalcMean

            #increment the counter by band...not needed to run the program
            CalcCount +=1
       

    #Write array of caclulated means to image file (ie. for a given band)
    #Notice this is indented to match the "band for loop"
            
    out_band.WriteArray(out_array)

    
    print "Created band: " + str(bnd)
    print "Total pixels Calculated: " + str(CalcCount)

#Clean up data files
#This may be necessary to see changes in output file
#if you need to redo
    
del outputimg
del img