High Dynamic Range Imaging

Name: Saaz Ramsubhag

Email:

Abstract

This project attempts to convert a series of multi-exposure images of the same scene to a displayable "HDR" image. This is done by recovering the camera's response function and subsequently generate the radiance map of the scene. We then consider different methods of tone mapping the high dynamic range radiance map to a displayable low dynamic range.

The camera response function and radiance map are recovered using the Debevec-Malik method, which allows us to solve for the camera response function and radiance of pixels as a overconstrained system of equations. Moreover, while Debevec & Malik outline the equations under the assumption that the exact exposure is known, we can loosen this constraint to only require an approximate exposure.

Lastly, we consider differnt methods of global tone mapping to display the HDR image without losing detail.

Code Libraries

The following basic code libraries are required to run this project: numpy, matplotlib and skimage.color.

• Numpy is heavily relied on for operations over multidimensional arrays as well as matrix operations.
• Matplotlib is used to display the images
• skimage.color is used to convert color images to grayscale

%matplotlib inline

import numpy as np
import numpy.linalg as la
import matplotlib
import matplotlib.image as image
import matplotlib.pyplot as plt
from skimage.color import rgb2gray

Auto Bracketed multi-exposure test images (grayscale)

def readImages(imgNames):
    imgs = []
    for n in imgNames:
        imgs.append(image.imread(n))

    return  np.array(imgs)

def grayscale(imgs):
    imgsgray = []
    for img in imgs:
        imgsgray.append(rgb2gray(img))
        imgsgray[-1] *= 255
    return np.round(np.array(imgsgray)).astype(int)

# HDR parameters
z_max = 255
z_min = 0
z_mid = (z_max + z_min) // 2
z_range = z_max - z_min + 1

# loading images and getting grayscale
cannonImgs = readImages(["images/cannon1.jpg","images/cannon2.jpg","images/cannon3.jpg"])
cannonTimes = np.array([1/350,1/180, 1/90])

imgs = cannonImgs
imgsgray = grayscale(imgs)
expTimes = cannonTimes
lnExpTimes = np.log(expTimes)
expNum = len(imgs)

best = image.imread("images/cannon.jpg")

fig = plt.figure(figsize = (15, 4))
plt.subplot(131)
plt.imshow(imgsgray[0], cmap = "gray")

plt.subplot(132)
plt.imshow(imgsgray[1], cmap = "gray")

plt.subplot(133)
plt.imshow(imgsgray[2], cmap = "gray")

<matplotlib.image.AxesImage at 0x2575fc851f0>

Recovering Hight Dynamic Range Radiance Map

The camera's radiometic response function, $f$, maps the exposure at a pixel's value to value of that pixel. Exposure is given by the product of the irradiance of that pixel and the exposure time the image was taken with. Hence $$z_{ij} = f(E_i t_j) $$ Where $z_{ij}$ is value of the pixel at location $i$ in the image with the $j$ th exposure time, $E_i$ is the irradiance at location $i$, and $j$ is the $j$ th exposure time. Instead of finding $f$, we will recover $g$ where $$ g(z_{ij}) = \ln f^{-1} = \ln E_i + \ln t_j$$

Response for assumed irradiance

Assuming the irradiance of every pixel is 1, $\ln E = 0$. Thus $g(z_{ij}) =\ln t_j$ and so we can plot the response for any pixel.

# random points in the images
sample_N = 4
markers = ['.-', '+-', 'x-', '*-']

sample = np.random.rand(sample_N, 2)
sample[:, 0] = sample[:, 0] * imgsgray[0].shape[0]
sample[:, 1] = sample[:, 1] * imgsgray[0].shape[1]
sample = np.round(sample).astype(int)

# plot response of pixels where each marker represents a pixel for different exposure times
fig, ax = plt.subplots()
for i in range(sample_N):
    ax.plot(imgsgray[:, sample[i, 0], sample[i, 1]], lnExpTimes, markers[i])
    
ax.set_xlabel('Pixel value (Zij)')
ax.set_ylabel('ln exposure(1 * tj)')
plt.title("g(z) for random points")
plt.show()

By plotting the response of each pixel over the different exposure times, we can see that each pixel gives a segment of the response curve. Thus by estimating the true irradiance values of the pixels, we can "shift" the curves to recover the true response curve

Estimating $g(z)$, $E_i$ and $\ln t_j$

Since we do not know $g(z_{ij})$ and $E_i$, we estimate both. Note that since $z$ refers to a pixel value, the range of values for $z$ is finite and discrete. For a given image, we assume that the minimum value, $z_{min} = 0$ and the max $z_{max} = 255$. Thus $g(z)$ has a finite range that we can solve for.

Additionally, we do not need to assume that the given exposure times are the exact exposure times of the images. Rather we can estimate $t_j$ such that it is close to the given exposure time $\hat t_j$.

To recover $g(z), E_i$ and $t_j$ such that they satisfy the equation for $g(z)$, the following least-squared error objective function is used $$O = \sum_{i=1}^N \sum_{j=1}^{P} \left[w(z_{ij})(g(z_{ij}) - \ln E_i - \ln t_j)\right]^2 + \lambda \sum_{y=z_{min+1}}^{z_{max-1}} [w(y)g''(y)]^2 + \eta \sum_{j=1}^P (\ln t_i - \ln \hat t_i)^2$$

Where $N$ and $P$ are the number of pixels and exposure times used respectively.

In addition to the least-squared error of $g$, each term of the sum is weighted by $w(z)$. Because of the anticipated shape of the response function, the data is expected to have a worse fit for $g(z)$ as $z$ nears $z_{min}$ or $z_{max}$. So the weight function emphasizes fitting at the middle of the curve $$ w(z) = \begin{cases} z - z_{min}, & z \leq \frac{1}{2}(z_{min} + z_{max}) \\ z_{max} - z, & z > \frac{1}{2}(z_{min} + z_{max}) \end{cases} $$

Additionally, $\lambda \sum_{y=z_{min+1}}^{z_{max-1}} [w(y)g''(y)]^2 $ is a smoothing term to enforce a zero second derivative where $$ g''(z) = g(z-1) - 2g(z) + g(z + 1)$$

The term $\eta \sum_{j=1}^P (\ln t_i - \ln \hat t_i)^2$ ensures that the estimated $\ln t_j$ is close to the given exposure time

Finally, the resulting $g$ and $E$ are "shift ambiguous" (i.e. their absolute values do not matter, only their values relative to each other). We remove this ambiguity with the condition that $g[z_{mid}] = 0$ where $z_{mid} = \frac{1}{2}(z_{max} + z_{min})$

# weight function
def Weight(z):
    if (z <= z_mid):
        return z - z_min
    return z_max - z

# store as array
w = np.fromfunction(np.vectorize(Weight), (z_max - z_min + 1,), dtype=int)

# returns the values of N randomly selected points from the images
#    returns a 1D array such that the first N elements are the values of points for the first image,
#    the second N elements are the values of the second image and so on
#    Does not return pixels that are fully saturated (0 or 255) for all exposures
def Z(imgs, N):
    P = len(imgs)
    points = np.empty((0,2), dtype=int)
    n = 0

    # discarding and resampling points at full saturation
    while n < N:
        p = np.random.rand(N-n, 2)
        p[:, 0] = p[:, 0] * imgs[0].shape[0]
        p[:, 1] = p[:, 1] * imgs[0].shape[1]
        p = np.floor(p).astype(int)

        mask = np.logical_and(imgs[:,p[:,0],p[:,1]] != np.full((P,1),255),
                              imgs[:,p[:,0],p[:,1]] != np.full((P,1),0))
        mask = np.logical_or.reduce(mask, 0)
        
        p = p[mask]
        n += len(p)
        
        points = np.concatenate((points, p), axis=0)
    
    z = imgs[:, points[:,0], points[:,1]]

    return z.flatten()

Minimum points needed

We don't need to use all pixels in an image to find $g$, but for a sufficently overdetermined system $$N(P-1) > (z_{max} - z_{min}) $$ Thus the minimum number of points needed is $$N_{min} = \left\lceil \frac{z_{max} - z_{min}}{P-1} \right\rceil $$

Matrix Form

To solve for $g$, $E_i$ and $t_j$, we can represent the problem as an overdetermined system of linear equations in matrix form Ax = b. Where we solve for $x = A^{-1}b$ using SVD

Where $x$ is the vector $\left[g(0),\ldots, g(255), \ln E_1, \ldots, \ln E_N, \ln t_1, \ldots, \ln t_P\right]^T$

We structure the matrix $A$ and vector $b$ as below

For the first $N*P$ rows, multiplying by $x$ gives the equation $w(z_{ij})(g(z_{ij}) - \ln E_i - \ln t_j) = 0$ for that point. For each row, the only value in Curve weights is $w(z)$ where $z$ is the value of that pixel and it is at the $zth$ column. The only value in RadWeights is $-w(z)$ at the column of the point it corresponds to. The only value in TimeWeights is $-w(z)$ at the column of the image that pixel corresponds to.

For the next $254$ rows, multiplying by $x$ gives the equation $\lambda w(z) \left( g(z-1) - 2g(z) + g(z + 1)\right) = 0$ for that value of $z$. The three values in each row of SmoothWeights are $\lambda w(z), -2\lambda w(z), \lambda w(z)$ where $z$ is the row and column of the middle value.

Multiplying the next $P$ rows gives the equation $\eta (\ln t_i - \ln \hat t_i) = 0$. Where its a diagonal of $\eta$

Finally, multiplying the last row gives $w(z_{mid})g(z_{mid}) = 0$. The only value is $w(z_{mid})$ at column $z_{mid}$

# solves for x where Ax = b using the SVD for A
def solveSVD(A,b):
    U, W, V_t = la.svd(A, full_matrices=False)
    W = np.diag(W)
    x = np.matmul(np.transpose(U), b)
    Wi = la.inv(W)
    x = np.matmul(Wi, x)
    x = np.matmul(np.transpose(V_t), x)
    return x

# estimate response function g, radiance, ln E and exposure time, ln t for the given pixels z
# from N points and P exposures
def EstimateResponse(z, lnt, N, P, lmda=1.0, eta=10.0):
    # CurveWeights (N*P, z_range):
    #  each row refers to the pixel index i+(j*N), the singular non zero entry of the
    #  row is w(z) at the column z where z is the value of the pixel
    x, y = np.ogrid[:N*P, :z_range]
    CurveWeights = np.where(z[x]==y, w[z[x]], 0)

    # RadWeights (N*P, N):
    #  each row refers to the pixel index i+(j*N), the singular non zero entry of the
    #  row is -w(z) at the column i where z is the value of the pixel
    x, y = np.ogrid[:N*P, :N]
    RadWeights = np.where(x%N==y, -1 * w[z[x]], 0)

    # TimeWeights (N*P, P)
    x, y = np.ogrid[:N*P, :P]
    TimeWeights = np.where(x//N==y, -1 * w[z[x]], 0)

    # SmoothWeights (z_range-2, z_range):
    #  tridiagonal matrix with the diagonal entries being lambda*w(z), -2*lambda*w(z), lambda*w(z)
    #  where z is the column of the middle entry
    d1 = np.arange(1,z_range-1)
    d1 = w[d1]
    d1 = np.diag(d1, k=2).astype(float)
    d1 = lmda * d1[0:z_range-2]

    d2 = np.arange(1,z_range)
    d2 = w[d2]
    d2 = np.diag(d2, k=1).astype(float)
    d2 = -2 * lmda * d2[0:z_range-2]

    d3 = np.arange(1,z_range+1)
    d3[-1] = 0
    d3 = w[d3]
    d3 = np.diag(d3, k=0).astype(float)
    d3 = lmda * d3[0:z_range-2]
    
    SmoothWeights = d2 + d1 + d3

    # Zeros1 (z_range-2, N+P):
    Zeros1 = np.zeros((z_range-2, N+P))

    x, y = np.ogrid[:P, :P]
    Times = np.where(x==y, eta, 0)

    # Zeros2 (P, z_range+N):
    Zeros2 = np.zeros((P, z_range+N))

    # MID (1, z_range + N + P): Add row to constrain g[z_mid] = 0
    Mid = np.zeros((1,z_range+N+P))
    Mid[0,z_mid] = w[z_mid] 

    # Joining submatrices
    A = np.concatenate((np.concatenate((CurveWeights, RadWeights, TimeWeights), axis=1), 
                        np.concatenate((SmoothWeights,Zeros1), axis=1),
                        np.concatenate((Zeros2,Times), axis=1),
                        Mid), axis=0, dtype=float)

    # b is a vector of exposure times, padded with zeros
    b = np.zeros(shape=(N*P+z_range-2+P+1))
    for j in range(P):
        b[N*P+z_range-2 + j] = eta*lnt[j]

    x = solveSVD(A, b)

    g = x[0:z_range]
    lnE = x[z_range:z_range+N]
    estlnt = x[z_range+N:]

    return g, lnE, estlnt

Plotting recovered E and g

# Estimate g and E
P = expNum
N = np.ceil((z_max - z_min) / (P-1)).astype(int)
z = Z(imgsgray, N)
g, lnE, estLnTimes= EstimateResponse(z=z, lnt=lnExpTimes, N=N, P=P, lmda=6, eta=20.0)

print('Estimated Exposures differences: {}'.format(lnExpTimes - estLnTimes))

# plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,4))

for i in range(N*P):
    ax1.plot(z[i], lnE[i%N] + estLnTimes[i//N], '.k')
    
ax1.set_xlabel('Pixel value (Zij)')
ax1.set_ylabel('log exposure(Ei * ti)')

ax2.plot(np.arange(z_min,z_max+1), g)   
ax2.set_xlabel('pixel value z')
ax2.set_ylabel('g(z)')
fig.suptitle('Response Curve from Recovered E, t and g')


plt.show()

Estimated Exposures differences: [-0.2427796  -0.01007908  0.25285868]

Color image

So far we have recovered the response curve for a grayscale image. For a color image, we recover the response curve for each of the RGB channels separately

def ColorZ(imgs, N):
    zRed   = Z(imgs[:,:,:, 0], N)
    zGreen = Z(imgs[:,:,:, 1], N)
    zBlue  = Z(imgs[:,:,:, 2], N)

    return np.array([zRed,zGreen,zBlue,])

def ColorResponse(zRGB, lnt, N, P, lmda=5.0, eta=10.0):
    gRed,   lnERed,   lnTRed   = EstimateResponse(z=zRGB[0], lnt=lnt, N=N, P=P, lmda=lmda, eta=eta)
    gGreen, lnEGreen, lnTGreen = EstimateResponse(z=zRGB[1], lnt=lnt, N=N, P=P, lmda=lmda, eta=eta)
    gBlue,  lnEBlue,  lnTBlue  = EstimateResponse(z=zRGB[2], lnt=lnt, N=N, P=P, lmda=lmda, eta=eta)

    gRGB = np.array([gRed,gGreen,gBlue])
    lnERGB = np.array([lnERed, lnEGreen, lnEBlue])
    lnTRGB = np.array([lnTRed,lnTGreen,lnTBlue])

    return gRGB, lnERGB, lnTRGB

P = expNum
N = np.ceil((z_max - z_min) / (P-1)).astype(int) 
zRGB = ColorZ(imgs, N)
gRGB, lnERGB, lnTRGB = ColorResponse(zRGB=zRGB, lnt=lnExpTimes, N=N, P=P, lmda=6, eta=20.0)

fig, axs = plt.subplots(1, 3, figsize=(12, 4))
for j in range(P):
    axs[0].plot(zRGB[0, j*N : (j+1)*N], lnERGB[0] + lnTRGB[0,j], '.r')
axs[0].set_title('Red')

for j in range(P):
    axs[1].plot(zRGB[1, j*N : (j+1)*N], lnERGB[1] + lnTRGB[1,j], '.g')
axs[1].set_title('Green')

for j in range(P):
    axs[2].plot(zRGB[2, j*N : (j+1)*N], lnERGB[2] + lnTRGB[2,j], '.b')
axs[2].set_title('Blue')


for ax in axs.flat:
    ax.set(xlabel='Pixel value (Zij)', ylabel='ln exposure(Ei * tj)')

for ax in axs.flat:
    ax.label_outer()

Generate Radiance Map

Now that we can recover the response curve of the camera, we can use it to generate a radiance map for each pixel in an image taken with it using $$ \ln E_i = g(z_{ij}) - \ln t_j $$ For robustness, we take a weighted average over the available exposures where the middle of the curve has the higher weight $$ \ln E_i = \frac{\sum_{j=1}^P w(z_{ij})(g(z_{ij}) - \ln t_j )}{\sum_{j=1}^P w(z_{ij})}$$

def RadianceMap(imgs, g, lnt):
    P = len(imgs)

    bottom = np.zeros(shape=imgs[0].shape)
    top = np.zeros(shape=imgs[0].shape)

    for p in range(P):
        bottom += w[imgs[p]]
        top += w[imgs[p]] * (g[imgs[p]] - lnt[p])

    # Edge case for when w(zij) == 0 for all j. Set value instead of averaging
    fullBright = np.argwhere(np.logical_and(bottom == 0, imgs[0] == z_max))
    fullDark = np.argwhere(np.logical_and(bottom == 0, imgs[0] == z_min))
    top[fullBright[:,0], fullBright[:,1]] = g[z_max] - np.min(lnt)
    top[fullDark[:,0], fullDark[:,1]] = g[z_min] - np.max(lnt)
    bottom[bottom == 0] = 1

    lnRadMap = np.divide(top,bottom)
    return np.exp(lnRadMap)

def ColorRadianceMap(imgs, gRGB, lnTRGB):
    radMapRGB = np.empty(shape=imgs[0].shape)
    for i in range(3):
        radMapRGB[:,:,i] = RadianceMap(imgs[:,:,:,i], gRGB[i], lnTRGB[i])

    return radMapRGB

Generate Displayable HDR image

The final radiance map tell us how bright the pixels in the image are relative to the others, giving us a High Dynamic Range image. To actually display the image, we must map the high dynamic range image to the low dynamic range of a display. The false-color image below shows the relative radiance values for the grayscale image.

radMap = RadianceMap(imgsgray, g, estLnTimes)
radMapRGB = ColorRadianceMap(imgs, gRGB, lnTRGB)
print('Minimum irradiance: {}'.format(np.min(radMap)))
print('Maximum irradiance: {}'.format(np.max(radMap)))

fig = plt.figure(figsize=(12,4))
plt.imshow(radMap, cmap='gray')
plt.title("Radiance map")
plt.show()

fig = plt.figure(figsize=(12,4))
plt.imshow(radMap, cmap='nipy_spectral')
plt.colorbar()
plt.title("False-color radiance map")
plt.show()

Minimum irradiance: 49.00938436250708
Maximum irradiance: 472.61889853024604

We can see that the range of values spans an order of magnitude and that most of the values lie at the bottom of the range. Thus simply linearly mapping the values to the low dynamic range of a display loses detail in the darker areas

Tone Mapping

Tone mapping is the process of mapping a radiance map to a low dynamic range for the purposes of displaying the image. As section 10.2.1 covers in the text book, global tone mapping is the simplest method. Global tone mapping is simply point processing.

One method of global tone mapping is taking the logarithm of the radiance map, and another is to apply gamma correction. For gamma correction, we use $$ I_{out} = \left(\frac{I_{in} - min(I_{in})}{max(I_{in}) - min(I_{in})}\right)^{\gamma}$$ Which maps intensity $I_{in}$ from the domain $[min(I_{in}), max(I_{in})]$ to the range $[0,1]$

def gammaCorrection(radMap, gamma):
    max_r = np.max(radMap)
    min_r = np.min(radMap)
    return ((radMap-min_r)/(max_r - min_r))** gamma

fig = plt.figure(figsize = (15, 4))
plt.subplot(131)
plt.imshow(gammaCorrection(radMap, gamma=0.2), cmap = "gray")
plt.title('gamma = 0.2')

plt.subplot(132)
plt.imshow(gammaCorrection(radMap, gamma=0.4), cmap = "gray")
plt.title('gamma = 0.4')

plt.subplot(133)
plt.imshow(gammaCorrection(radMap, gamma=0.8), cmap = "gray")
plt.title('gamma = 0.8')

Text(0.5, 1.0, 'gamma = 0.8')

We can see that as $\gamma$ decrease, we can see more detail in the darker areas of the image, but we lose detail in the brighter regions and the image appears more "washed out"

fig = plt.figure(figsize = (15, 13))
plt.subplot(231)
plt.imshow(imgsgray[1], cmap = "gray")
plt.title('Middle Exposure')

plt.subplot(232)
plt.imshow(rgb2gray(best), cmap = "gray")
plt.title('Example HDR')

plt.subplot(131)
plt.imshow(np.log(radMap), cmap = "gray")
plt.title('Log Tone Mapped')

plt.subplot(132)
plt.imshow(gammaCorrection(radMap, gamma=0.7), cmap = "gray")
plt.title('Gamma Correction')

Text(0.5, 1.0, 'Gamma Correction')

We can see that gamma correction is the better of the two methods while preserving detail throughout the image, but it loses contrast, making the image appear washed out compared to the original middle exposure or the example HDR image

Color Tone Mapping

From the textbook, we can use gamma correction to tone map color images in two ways: applying gamma correction to each of the RGB channels or by extracting the luminance of the radiance map, applying gamma correction to that and setting that to be the new luminance of the image.

The luminance, Y, of an RGB image is given by $$Y = 0.17697R + 0.81240G + 0.01063B$$ where the image is in the domain [0,1]

To set the luminance of the original image with the new luminance $L_{out}$, we use $$C_{out} = \left( \frac{C_{in}}{L_{in}} \right)^s L_{out}$$ where $C_{in} = (R,G,B)$ and $L_{in}$ are the original image and luminance and $C_{out}$ and $L_{out}$ are the values after correction

# Linearly maps each channel of the img to the range [0,1] separately
def ColorLinMap(radMapRGB):
    map = radMapRGB.copy()
    for i in range(3):
        map[:,:,i] = np.interp(map[:,:,i], [np.min(map[:,:,i]), np.max(map[:,:,i])], [0, 1])
    return map

# Linearly maps the img from the domain [0,1] to the range [0,255] as integers
def IntegerColor(img):
    map = img.copy()
    for i in range(3):
        map[:,:,i] = np.interp(map[:,:,i], [0, 1], [0, 255])
    return np.round(map).astype(int)

# Performs gamma correction on each channel separately to the range [0,1]
def BasicColorGamma(radMapRGB, gamma):
    map = radMapRGB.copy()
    for i in range(3):
        map[:,:,i] = gammaCorrection(map[:,:,i], gamma)
    return map

# Extract the luminance from an rgb img in the domain [0,1]
def ExtractLum(img):
    return np.dot(img, np.array([0.17697, 0.81240, 0.01063]))

# Set the luminance for an rgb img in the domain [0,1]
def SetLum(img, newLum, s=0.6):
    oldLum = ExtractLum(img)
    return ((np.divide(img, np.expand_dims(oldLum,axis=2)))**s) * np.expand_dims(newLum,axis=2)

# Performs gamma correction on a radiance map to an img in the range [0,1]
def ColorGamma(radMapRGB, gamma):
    normImg = ColorLinMap(radMapRGB)
    lum = ExtractLum(normImg)
    lum = lum ** gamma
    return SetLum(normImg, lum)

fig = plt.figure(figsize = (15, 13))
plt.subplot(231)
plt.imshow(imgs[1])
plt.title('Middle Exposure')

plt.subplot(232)
plt.imshow(best)
plt.title('easyHDR Sample')

plt.subplot(131)
plt.imshow(IntegerColor(BasicColorGamma(radMapRGB, 0.65)))
plt.title('Per Channel Gamma Correction')

plt.subplot(132)
plt.imshow(IntegerColor(ColorGamma(radMapRGB, 0.65)))
plt.title('Luminance Gamma Correction')

Text(0.5, 1.0, 'Luminance Gamma Correction')

Again, we can see that the tone mapping preserves the details of both dark and bright areas of the image but it appears much more "washed out" compared to the original middle exposure or example HDR image.

Additionally, while the Luminance gamma correction appears slightly more saturated than the basic gamma correction, the difference is marginal.

More examples

def GenerateMappedHDR(imgs, exposureTimes, lmda=5.0, eta=50.0, gamma=0.6):
    P = expNum
    N = 2 * np.ceil((z_max - z_min) / (P-1)).astype(int)

    zRGB = ColorZ(imgs, N)
    gRGB, lnERGB, lnTRGB = ColorResponse(zRGB=zRGB, lnt=np.log(exposureTimes), N=N, P=P, lmda=lmda, eta=eta)
    radMapRGB = ColorRadianceMap(imgs, gRGB, lnTRGB)

    return IntegerColor(ColorGamma(radMapRGB, gamma))

blimpImgs = readImages(["images/blimp1.jpg","images/blimp2.jpg","images/blimp3.jpg"])
blimpBest = image.imread("images/blimp.jpg")
blimpTimes = np.array([1/350,1/250, 1/125])

sunsetImgs = readImages(["images/wadi-rum-sunset1.jpg","images/wadi-rum-sunset2.jpg","images/wadi-rum-sunset3.jpg"])
sunsetBest = image.imread("images/wadi-rum-sunset.jpg")
sunsetTimes = np.array([1/500,1/250, 1/60])

capeImgs = readImages(["images/cap-de-formentor1.jpg","images/cap-de-formentor2.jpg","images/cap-de-formentor3.jpg"])
capeBest = image.imread("images/cap-de-formentor.jpg")
capeTimes = np.array([1/4000,1/1500, 1/350])

fig = plt.figure(figsize = (15, 10))

plt.subplot(331)
plt.imshow(capeImgs[1])
plt.title('Middle Exposure')

plt.subplot(332)
plt.imshow(capeBest)
plt.title('easyHDR Sample')

plt.subplot(333)
plt.imshow(GenerateMappedHDR(capeImgs, capeTimes, lmda=6, gamma=0.35))
plt.title('Gamma Tone Mapping')

plt.subplot(334)
plt.imshow(blimpImgs[1])
plt.title('Middle Exposure')

plt.subplot(335)
plt.imshow(blimpBest)
plt.title('easyHDR Sample')

plt.subplot(336)
plt.imshow(GenerateMappedHDR(blimpImgs, blimpTimes, lmda=10, gamma=0.6))
plt.title('Gamma Tone Mapping')

plt.subplot(337)
plt.imshow(sunsetImgs[1])
plt.title('Middle Exposure')

plt.subplot(338)
plt.imshow(sunsetBest)
plt.title('easyHDR Sample')

plt.subplot(339)
plt.imshow(GenerateMappedHDR(sunsetImgs, sunsetTimes, lmda=10, gamma=0.5))
plt.title('Gamma Tone Mapping')

Text(0.5, 1.0, 'Gamma Tone Mapping')

fig = plt.figure(figsize = (15, 4))
plt.subplot(131)
plt.imshow(GenerateMappedHDR(capeImgs, capeTimes, lmda=6, gamma=0.35))

plt.subplot(132)
plt.imshow(GenerateMappedHDR(capeImgs, capeTimes, lmda=6, gamma=0.35))

plt.subplot(133)
plt.imshow(GenerateMappedHDR(capeImgs, capeTimes, lmda=6, gamma=0.35))

<matplotlib.image.AxesImage at 0x2577033bbf0>

Note that my implementation of tone mapping is inconsistent for the same input parameters where the dominant hue of the image varies greatly. This could be due to a bug in the implementation where the "strength" of one channel compared to the others randomly varies with small variations in the estamated response curve

Conclusions

In applying the Debevec-Malik method, we were able to recover smooth camera response curves from the multi-exposure images and the resulting radiance maps had relatively high dynamic ranges. Moreover, we were able to drop the requirement of knowing the exact exposure time of the images and replace it with an approximate exposure time. For color images, we ended up with 3 different sets of exposure times, generated from each of the channels separately, which could have negatively affected the results. One way to ensure that exposure times are the same across the 3 channels would be to solve for their response curves together in the same least squares problem.

Looking at the results of the global tone mapping, we see that high-level detail is preserved across dark and light areas but the fine details are not, making the result washed out. One way to avoid this could be to use local tone mapping to take into account how the radiance of a pixel differs from that of its neighbours. Additionally, the gamma correction of color images appears to be unstable. One possible solution is change how I normalize the color channels repective to each other

While there are some issues, I was able to perform high dynamic range imaging over multiple exposures.