From 77b05cc86cf333181edffeb3917a849cca302c72 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 21 Oct 2022 15:10:14 -0400 Subject: [PATCH] Wrote an override of the matplotlib LightSource class to remove the intensity scaling step which was causing the movies of the surface evolution to flicker. --- python/ctem/ctem/util.py | 69 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 63 insertions(+), 6 deletions(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 488280dc..c750e774 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -5,6 +5,62 @@ import matplotlib.cm as cm import matplotlib.pyplot as plt + +class CTEMLightSource(LightSource): + """ + Override one function in LightSource to prevent the contrast from being rescaled. + """ + def shade_normals(self, normals, fraction=1.): + """ + Calculate the illumination intensity for the normal vectors of a + surface using the defined azimuth and elevation for the light source. + + Imagine an artificial sun placed at infinity in some azimuth and + elevation position illuminating our surface. The parts of the surface + that slope toward the sun should brighten while those sides facing away + should become darker. + + Changes by David Minton: The matplotlib version of this rescales the intensity + in a way that causes the brightness level of the images to change between + simulation outputs. This makes movies of the surface evolution appear to flicker. + + Parameters + ---------- + fraction : number, optional + Increases or decreases the contrast of the hillshade. Values + greater than one will cause intermediate values to move closer to + full illumination or shadow (and clipping any values that move + beyond 0 or 1). Note that this is not visually or mathematically + the same as vertical exaggeration. + + Returns + ------- + ndarray + A 2D array of illumination values between 0-1, where 0 is + completely in shadow and 1 is completely illuminated. + """ + + intensity = normals.dot(self.direction) + + # Apply contrast stretch + imin, imax = intensity.min(), intensity.max() + intensity *= fraction + + # Rescale to 0-1, keeping range before contrast stretch + # If constant slope, keep relative scaling (i.e. flat should be 0.5, + # fully occluded 0, etc.) + # if (imax - imin) > 1e-6: + # # Strictly speaking, this is incorrect. Negative values should be + # # clipped to 0 because they're fully occluded. However, rescaling + # # in this manner is consistent with the previous implementation and + # # visually appears better than a "hard" clip. + # intensity -= imin + # intensity /= (imax - imin) + intensity = np.clip(intensity, 0, 1) + + return intensity + + # Set pixel scaling common for image writing, at 1 pixel/ array element dpi = 72.0 @@ -125,16 +181,17 @@ def create_rplot(user, odist, pdist, tdist, ph1): return + def image_dem(user, DEM): dpi = 300.0 # 72.0 pix = user['pix'] gridsize = user['gridsize'] ve = 1.0 azimuth = 300.0 # user['azimuth'] - solar_angle = 15.0 # user['solar_angle'] + solar_angle = 20.0 # user['solar_angle'] - ls = LightSource(azdeg=azimuth, altdeg=solar_angle) - dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix) + ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle) + dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix, fraction=1) # Generate image to put into an array height = gridsize / dpi @@ -180,11 +237,11 @@ def image_shaded_relief(user, DEM): pix = user['pix'] gridsize = user['gridsize'] ve = 1.0 - mode = 'hsv' + mode = 'overlay' azimuth = 300.0 # user['azimuth'] - solar_angle = 15.0 # user['solar_angle'] + solar_angle = 20.0 # user['solar_angle'] - ls = LightSource(azdeg=azimuth, altdeg=solar_angle) + ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle) cmap = cm.cividis # If min and max appear to be reversed, then fix them