例如设计一个3 x 3的滑动窗口,写算法执行就有两种方式:
1.pixel by piexl每个进行逐像素运算,效率太低,速度慢
2.使用 slice切片形式循环,效率高,速度快
两个作业就是分别用pixel和slice方式完成高通滤波操作进行对比
1.Assignment 6a
- Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img
- The output data type will be Float
- Use pixel notation (that’s why you’re doing it on smallaster.img instead of aster.img)
1)自己写的代码,对输入图像的三个通道进行高通滤波,输出图像有三个图像
# Assignment 6a:Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img
# writen by myself
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys
import time
# record start time
startTime = time.time()
# set the directory
os.chdir('E:/data/GDAL/ospy_data6')
# register all driver
gdal.AllRegister()
# open file
ds = gdal.Open('smallaster.img', GA_ReadOnly)
if ds is None:
print('Cannot open smallaster.img')
sys.exit(1)
# get the size of smallaster.img
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
# create a data array to store result
outdata = np.zeros((bands, rows, cols))
# filt = np.full((3, 3), 0.111).astype(np.float32)# low pass
filt = np.array([[-0.7, -1.0, -0.7], [-1.0, 6.8, -1.0], [-0.7, -1.0, -0.7]]).astype(np.float32) # high pass
# get current driver
driver = ds.GetDriver()
# create output file
outds = driver.Create('highpass_smallaster.img', cols, rows, bands, GDT_Float32)
if outds is None:
print('Cannot open highpass_smallaster.img')
sys.exit(1)
# cycle bands to calculate
for i in range(bands):
# get each band
band = ds.GetRasterBand(i + 1)
# get each band data
data = band.ReadAsArray(0, 0, cols, rows).astype(np.float32)
# cycle row direction
for j in range(1, rows - 1):
# cycle col direction
for k in range(1, cols - 1):
# high pass filter calculate
outdata[i, j, k] = (
data[j - 1, k - 1] * filt[0, 0] + data[j - 1, k] * filt[0, 1] + data[j - 1, k + 1] * filt[0, 2] +
data[j, k - 1] * filt[1, 0] + data[j, k] * filt[1, 1] + data[j, k + 1] * filt[1, 2] +
data[j + 1, k - 1] * filt[2, 0] + data[j + 1, k] * filt[2, 1] + data[j + 1, k + 1] * filt[2, 2])
# get current band to store result
outband = outds.GetRasterBand(i + 1)
outband.WriteArray(outdata[i, :, :], 0, 0)
# write to disk
outband.FlushCache()
stats = outband.GetStatistics(0, 1)
# set the Geotransform and projection
outds.SetGeoTransform(ds.GetGeoTransform())
outds.SetProjection(ds.GetProjection())
del ds
del outds
# record the end time
endTime = time.time()
# print the time for run last
print('Last ', endTime - startTime, ' seconds')
程序运行时间:
Last 104.1204240322113 seconds
2)参考答案给的代码,对输入图像的第一个通道进行高通滤波,输出图像只有一个通道
# homework 6a:high-pass filter using pixel notation
# reference answer
import os, sys, numpy, gdal, time
from gdalconst import *
start = time.time()
os.chdir('E:/data/GDAL/ospy_data6')
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open('smallaster.img', GA_ReadOnly)
if inDs is None:
print('Could not open smallaster.img')
sys.exit(1)
# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize
# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(numpy.float)
# do the calculation
outData = numpy.zeros((rows, cols), numpy.float)
for i in range(1, rows - 1):
for j in range(1, cols - 1):
outData[i, j] = ((-0.7 * inData[i - 1, j - 1]) + (-1.0 * inData[i - 1, j]) + (-0.7 * inData[i - 1, j + 1]) +
(-1.0 * inData[i, j - 1]) + (6.8 * inData[i, j]) + (-1.0 * inData[i, j + 1]) +
(-0.7 * inData[i + 1, j - 1]) + (-1.0 * inData[i + 1, j]) + (-0.7 * inData[i + 1, j + 1]))
# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('highpass1.img', cols, rows, 1, GDT_Float32)
if outDs is None:
print('Could not create highpass1.img')
sys.exit(1)
outBand = outDs.GetRasterBand(1)
# write the output data
outBand.WriteArray(outData, 0, 0)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
inDs = None
outDs = None
print(time.time() - start, 'seconds')
程序运行时间:
34.041293144226074 seconds
2.Assignment 6b
- Use a 3x3 high pass filter to detect edges in band 1 of aster.img (good idea to test on smallaster.img first)
- The output data type will be Float
- Use slice notation
1)自己写的代码,输出图像有三个通道
# Assignment 6b:Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img
# writen by myself
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys
import time
# record start time
startTime = time.time()
# set the directory
os.chdir('E:/data/GDAL/ospy_data6')
# register all driver
gdal.AllRegister()
# open file
ds = gdal.Open('smallaster.img', GA_ReadOnly)
if ds is None:
print('Cannot open smallaster.img')
sys.exit(1)
# get the size of smallaster.img
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
# create a data array to store result
outdata = np.zeros((bands, rows, cols))
# filt = np.full((3, 3), 0.111).astype(np.float32)# low pass
filt = np.array([[-0.7, -1.0, -0.7], [-1.0, 6.8, -1.0], [-0.7, -1.0, -0.7]]).astype(np.float32) # high pass
# get current driver
driver = ds.GetDriver()
# create output file
outds = driver.Create('highpass_smallaster_slice.img', cols, rows, bands, GDT_Float32)
if outds is None:
print('Cannot open highpass_smallaster_slice.img')
sys.exit(1)
# cycle bands to calculate
for i in range(bands):
# get each band
band = ds.GetRasterBand(i + 1)
# get each band data
data = band.ReadAsArray(0, 0, cols, rows).astype(np.float32)
outdata[i, 1:rows - 1, 1:cols - 1] = (
data[0:rows - 2, 0:cols - 2] * filt[0, 0] + data[0:rows - 2, 1:cols - 1] * filt[0, 1] +
data[0:rows - 2, 2:cols] * filt[0, 2] + data[1:rows - 1, 0:cols - 2] * filt[1, 0] +
data[1:rows - 1, 1:cols - 1] * filt[1, 1] +
data[1:rows - 1, 2:cols] * filt[1, 2] + data[2:rows, 0:cols - 2] * filt[2, 0] +
data[2:rows, 1:cols - 1] * filt[2, 1] + data[2:rows, 2:cols] * filt[2, 2])
# get current band to store result
outband = outds.GetRasterBand(i + 1)
outband.WriteArray(outdata[i, :, :], 0, 0)
# write to disk
outband.FlushCache()
stats = outband.GetStatistics(0, 1)
# set the Geotransform and projection
outds.SetGeoTransform(ds.GetGeoTransform())
outds.SetProjection(ds.GetProjection())
del ds
del outds
# record the end time
endTime = time.time()
# print the time for run last
print('Last ', endTime - startTime, ' seconds')
# last 3s
程序运行时间,可以看到速度由104s提升到了3s
Last 3.429971933364868 seconds
2)参考答案给的代码,输出图形只有一个通道
# homework 6b:high-pass filter using slice notation
# reference answer
import os, sys, numpy, gdal, time
from gdalconst import *
start = time.time()
os.chdir('E:/data/GDAL/ospy_data6')
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open('aster.img', GA_ReadOnly)
if inDs is None:
print('Could not open aster.img')
sys.exit(1)
# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize
# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(numpy.float)
# do the calculation
outData = numpy.zeros((rows, cols), numpy.float)
outData[1:rows - 1, 1:cols - 1] = (
(-0.7 * inData[0:rows - 2, 0:cols - 2]) +
(-1.0 * inData[0:rows - 2, 1:cols - 1]) + (-0.7 * inData[0:rows - 2, 2:cols]) +
(-1.0 * inData[1:rows - 1, 0:cols - 2]) + (6.8 * inData[1:rows - 1, 1:cols - 1]) +
(-1.0 * inData[1:rows - 1, 2:cols]) + (-0.7 * inData[2:rows, 0:cols - 2]) +
(-1.0 * inData[2:rows, 1:cols - 1]) + (-0.7 * inData[2:rows, 2:cols]))
# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('highpass3.img', cols, rows, 1, GDT_Float32)
if outDs is None:
print('Could not create highpass3.img')
sys.exit(1)
outBand = outDs.GetRasterBand(1)
# write the output data
outBand.WriteArray(outData, 0, 0)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
inDs = None
outDs = None
print(time.time() - start, 'seconds')
数据下载:https://download.csdn.net/download/qq_37935516/10798222