[Image-SIG] Two PIL based Hill Shading Implementations, using Python & Python/C.

John Barratt jb at langarson.com.au
Thu Feb 8 12:37:47 CET 2007

Hi Will,

Will Henney wrote:
> How about this for a numpy implementation? Runs nearly as fast as the C version
> on my mac (1.58 vs 1.16 secs). 
This is great!  I get about 0.34 vs 0.55.  The numpy implementation is 
darn good for a python only implementation.  There is a subtle 
difference in the output though, this version seems to create a what 
appears as a slightly 'noisier' result, the python and C versions I did 
were identical to the eye at least at 1:1.  Not sure exactly why this 
would be perhaps a value is clipping somewhere?  Could be also there is 
a consistent error in my two implementations.

I upgraded my numpy to the latest version and now it works with the 
latest PIL, so I will have to revisit the other simple pixel setting 
test that didn't work at all for me before.

I have also looked at the performance python article 
(http://www.scipy.org/PerformancePython) that was forwarded through to 
me (thanks Sebastian!) which looks like it opens up yet more convenient 
performance options with the weave modules, particularly now that I have 
numpy working with PIL.

Lots to explore still it seems!

Thanks for the example,



> def hillShadeNumPy(filenameIn, filenameOut, scale=1.0, azdeg=315.0, altdeg=45.0):
>     ''' Create a hill shade version of the given image using numpy '''
>     from numpy import sin, cos, hypot, arctan, arctan2, pi
>     dScale = (1.0 * scale)
>     # convert alt, az to radians
>     az = azdeg*pi/180.0
>     alt = altdeg*pi/180.0
>     (img, data, imgS, dataS) = initialiseImages(filenameIn)
>     # get the image data as a numpy array of floats
>     adata = N.asarray(img).astype('float')
>     # gradient in x and y directions
>     dx, dy = N.gradient(adata/dScale)
>     slope = 0.5*pi - arctan(hypot(dx, dy))
>     aspect = arctan2(dx, dy)
>     c = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(-az - aspect - 0.5*pi)
>     c = N.where(c > 0.0, c*255.0, 0.0)
>     imgS = Image.fromarray(c.astype('uint8'))
>     imgS.save(filenameOut,'PNG')

John Barratt - www.langarson.com.au
          Python, Zope, GIS, Weather

More information about the Image-SIG mailing list