''' a) Perform bilinear interpolation of the missing pixels in each

color channel
R, G and B to reconstruct a full color image from RawImage1 with the
CFA
provided in ‘[Link]’ '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Load [Link] and [Link]
# ---------------------------
raw_data = [Link]('Copy of [Link]') # Load raw Bayer
data
mask_data = [Link]('Copy of [Link]') # Load
CFA mask

raw_image = raw_data['RawImage1'] # Bayer raw image

cfa_mask = mask_data['bayer1'] # CFA mask: 1=R, 2=G, 3=B

rows, cols = raw_image.shape

# Coordinate grids
x, y = [Link]([Link](cols), [Link](rows))

# ---------------------------
# Interpolation function
# ---------------------------
def interpolate_channel(mask_value):
""" Interpolate missing values for a single channel using bilinear
interpolation """
points = np.column_stack((x[cfa_mask==mask_value],
y[cfa_mask==mask_value]))
values = raw_image[cfa_mask==mask_value]
full_channel = griddata(points, values, (x, y), method='linear')
# bilinear interpolation
return full_channel

# ---------------------------
# Interpolate R, G, B channels
# ---------------------------
R = interpolate_channel(1)
G = interpolate_channel(2)
B = interpolate_channel(3)

# ---------------------------
# Combine channels and normalize
# ---------------------------
RGB = [Link]([R, G, B], axis=-1)
RGB = RGB / [Link](RGB) # Normalize for display
RGB = np.nan_to_num(RGB) # Replace NaNs with 0

# ---------------------------
# Display raw and reconstructed images side by side
# ---------------------------
[Link](figsize=(10,5))

# Raw Bayer Image

[Link](1, 2, 1)
[Link](raw_image / [Link](raw_image), cmap='gray') # Normalize
raw for display
[Link]('Raw Bayer Image')
[Link]('off')

# Reconstructed RGB Image

[Link](1, 2, 2)
[Link](RGB)
[Link]('Reconstructed RGB (Bilinear Interpolation)')
[Link]('off')

plt.tight_layout()
[Link]()

# ---------------------------
# Save reconstructed image as PNG
# ---------------------------
[Link]("Reconstructed_RGB_a.png", RGB)
print("Reconstructed RGB image saved as 'Reconstructed_RGB_a.png'")

Reconstructed RGB image saved as 'Reconstructed_RGB_a.png'

''' b) Perform bicubic interpolation of the missing pixels in each

color channel
R, G and B to reconstruct a full color image from RawImage1 with the
CFA
provided in ‘[Link]’. '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Load [Link] and [Link]
# ---------------------------
raw_data = [Link]('Copy of [Link]') # Load raw Bayer
data
mask_data = [Link]('Copy of [Link]') # Load CFA mask

raw_image = raw_data['RawImage1'] # Bayer raw image

cfa_mask = mask_data['bayer1'] # CFA mask: 1=R, 2=G, 3=B

rows, cols = raw_image.shape

x, y = [Link]([Link](cols), [Link](rows)) # Coordinate
grids

# ---------------------------
# Interpolation function for bicubic
# ---------------------------
def interpolate_channel(mask_value):
""" Interpolate missing values for a single channel using bicubic
interpolation """
points = np.column_stack((x[cfa_mask == mask_value], y[cfa_mask ==
mask_value]))
values = raw_image[cfa_mask == mask_value]
full_channel = griddata(points, values, (x, y), method='cubic') #
bicubic interpolation
return full_channel

# ---------------------------
# Interpolate R, G, B channels using bicubic
# ---------------------------
R = interpolate_channel(1)
G = interpolate_channel(2)
B = interpolate_channel(3)

# ---------------------------
# Combine channels and normalize
# ---------------------------
RGB_bicubic = [Link]([R, G, B], axis=-1)
RGB_bicubic = RGB_bicubic / [Link](RGB_bicubic) # Normalize for
display
RGB_bicubic = np.nan_to_num(RGB_bicubic) # Replace NaNs
with 0

# ---------------------------
# Display raw and bicubic reconstructed images
# ---------------------------
[Link](figsize=(10,5))

# Raw Bayer Image

[Link](1, 2, 1)
[Link](raw_image / [Link](raw_image), cmap='gray') # Normalize
raw for display
[Link]('Raw Bayer Image')
[Link]('off')

# Bicubic Reconstructed RGB Image

[Link](1, 2, 2)
[Link](RGB_bicubic)
[Link]('Reconstructed RGB (Bicubic Interpolation)')
[Link]('off')

plt.tight_layout()
[Link]()

# Clip values to [0,1] to avoid errors

RGB_bicubic_safe = [Link](RGB_bicubic, 0, 1)

# Save the image

[Link]("Reconstructed_RGB_Bicubic_1_b.png", RGB_bicubic_safe)
print("Bicubic reconstructed RGB image saved as
'Reconstructed_RGB_Bicubic_1_b.png'")

WARNING:[Link]:Clipping input data to the valid range for

imshow with RGB data ([0..1] for floats or [0..255] for integers). Got
range [-0.03846349901167554..1.0].

Bicubic reconstructed RGB image saved as

'Reconstructed_RGB_Bicubic_1_b.png'
''' Compare the reconstructed full colour image with
that obtained using bilinear interpolation. '''

import numpy as np
import [Link] as plt
from [Link] import loadmat
from [Link] import griddata
from [Link] import sobel

# ---------- Load Data ----------

raw_data = loadmat('Copy of [Link]')
bayer_data = loadmat('Copy of [Link]')

raw_image = raw_data['RawImage1']
bayer_mask = bayer_data['bayer1']

H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

# ---------- Separate channels ----------

R = [Link](bayer_mask == 1, raw_image, [Link])
G = [Link](bayer_mask == 2, raw_image, [Link])
B = [Link](bayer_mask == 3, raw_image, [Link])

def interpolate_channel(channel, method):

known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method=method).reshape(H, W)

# ---------- Bilinear & Bicubic ----------

R_bilinear = interpolate_channel(R, 'linear')
G_bilinear = interpolate_channel(G, 'linear')
B_bilinear = interpolate_channel(B, 'linear')

R_bicubic = interpolate_channel(R, 'cubic')

G_bicubic = interpolate_channel(G, 'cubic')
B_bicubic = interpolate_channel(B, 'cubic')

bilinear_img = [Link]((R_bilinear, G_bilinear, B_bilinear)).clip(0,

255).astype(np.uint8)
bicubic_img = [Link]((R_bicubic, G_bicubic, B_bicubic)).clip(0,
255).astype(np.uint8)

# ---------- Sharpness Metric using Sobel ----------

def sharpness_metric(img):
gray = [Link](img, axis=2) # convert to grayscale
return [Link](sobel(gray)) # mean edge strength

sharp_bilinear = sharpness_metric(bilinear_img)
sharp_bicubic = sharpness_metric(bicubic_img)

print(f"Sharpness (Bilinear): {sharp_bilinear:.4f}")

print(f"Sharpness (Bicubic): {sharp_bicubic:.4f}")

if sharp_bicubic > sharp_bilinear:

print("Bicubic is sharper, likely better for edges.")
else:
print("Bilinear is smoother, may be better if noise is an issue.")

# ---------- Display ----------

[Link](figsize=(10,5))
[Link](1,2,1); [Link](bilinear_img); [Link]("Bilinear
Interpolation"); [Link]("off")
[Link](1,2,2); [Link](bicubic_img); [Link]("Bicubic
Interpolation"); [Link]("off")
plt.tight_layout()
[Link]()

# ---------- Save Images ----------

combined = [Link]((bilinear_img, bicubic_img))
[Link]("Reconstructed_Comparison.png", combined)

print("Saved 'Reconstructed_Comparison.png' (Bilinear | Bicubic)")

/tmp/[Link]: RuntimeWarning: invalid value

encountered in cast
bilinear_img = [Link]((R_bilinear, G_bilinear,
B_bilinear)).clip(0, 255).astype(np.uint8)
/tmp/[Link]: RuntimeWarning: invalid value
encountered in cast
bicubic_img = [Link]((R_bicubic, G_bicubic, B_bicubic)).clip(0,
255).astype(np.uint8)

Sharpness (Bilinear): 6.1259

Sharpness (Bicubic): 6.6290
Bicubic is sharper, likely better for edges.
Saved 'Reconstructed_Comparison.png' (Bilinear | Bicubic)

Aspect Bilinear Interpolation Bicubic Interpolation

Computation Simple, fast, uses linear kernel More complex, uses cubic kernel →
slower
Image Produces smoother images with less Preserves sharper edges, but may
Smoothness noise amplification cause ringing artifacts
Sharpness Lower sharpness → fewer details but Higher sharpness → retains more
Metric fewer artifacts details but amplifies noise
Noise Handles noisy images better due to May amplify noise if the image is
Sensitivity smoothing already noisy
Edge Blurs edges slightly Retains edge contrast better
Preservation
Use Case When noise reduction and When details and edges are
Preference smoothness are preferred important
''' c) Do demosaicing for [Link] using pattern ‘rggb’ with
matlab
built-in function ‘demosaic’ (or equivalently, [Link] with the
appropriate
Bayer conversion code in Python) and compare with previous two
reconstructed
images '''

import numpy as np
import [Link] as sio
import [Link] as plt
import cv2
from [Link] import griddata
from [Link] import sobel

# ---------- 1. Load raw image and CFA ----------

raw_data = [Link]('Copy of [Link]')
bayer_data = [Link]('Copy of [Link]')

raw_image = raw_data['RawImage1'].astype(np.uint16)
bayer_mask = bayer_data['bayer1']

H, W = raw_image.shape

# ---------- 2. Prepare interpolation ----------

X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

def interpolate_channel(channel, method):

known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method=method).reshape(H, W)

R = [Link](bayer_mask==1, raw_image, [Link])

G = [Link](bayer_mask==2, raw_image, [Link])
B = [Link](bayer_mask==3, raw_image, [Link])

# ---------- 3. Bilinear & Bicubic ----------

# Bilinear
R_lin = interpolate_channel(R, 'linear')
G_lin = interpolate_channel(G, 'linear')
B_lin = interpolate_channel(B, 'linear')
bilinear_img = [Link]((R_lin, G_lin,
B_lin)).clip(0,255).astype(np.uint8)

# Bicubic
R_cub = interpolate_channel(R, 'cubic')
G_cub = interpolate_channel(G, 'cubic')
B_cub = interpolate_channel(B, 'cubic')
bicubic_img = [Link]((R_cub, G_cub,
B_cub)).clip(0,255).astype(np.uint8)

# ---------- 4. Demosaicing using OpenCV ----------

raw_uint8 = raw_image.astype(np.uint8)
demosaic_img = [Link](raw_uint8, cv2.COLOR_BayerRG2BGR) # RGGB
pattern
demosaic_img = [Link](demosaic_img, cv2.COLOR_BGR2RGB)

# ---------- 5. Metrics ----------

def sharpness_metric(img):
gray = [Link](img, cv2.COLOR_RGB2GRAY)
return [Link](sobel(gray))

def noise_estimation(img):
gray = [Link](img, cv2.COLOR_RGB2GRAY)
return [Link](gray[100:150, 100:150]) # flat region

methods = {
"Bilinear": bilinear_img,
"Bicubic": bicubic_img,
"OpenCV Demosaic": demosaic_img
}

metrics = {}
for name, img in [Link]():
metrics[name] = {
"Sharpness": sharpness_metric(img),
"Noise": noise_estimation(img)
}

print("Image Quality Metrics (Higher Sharpness → Better Edges, Lower

Noise → Cleaner Image):")
for name, vals in [Link]():
print(f"{name:15s} | Sharpness: {vals['Sharpness']:.3f} | Noise:
{vals['Noise']:.3f}")

# ---------- 6. Display images ----------

[Link](figsize=(15,5))
[Link](1,3,1); [Link](bilinear_img); [Link]("Bilinear
Interpolation"); [Link]('off')
[Link](1,3,2); [Link](bicubic_img); [Link]("Bicubic
Interpolation"); [Link]('off')
[Link](1,3,3); [Link](demosaic_img); [Link]("OpenCV
Demosaic (RGGB)"); [Link]('off')
plt.tight_layout()
[Link]()

# ---------- 7. Save all images in one PNG ----------

fig, axes = [Link](1,3, figsize=(15,5))
for ax, img, title in zip(axes, [bilinear_img, bicubic_img,
demosaic_img],
["Bilinear", "Bicubic", "OpenCV Demosaic"]):
[Link](img)
ax.set_title(title)
[Link]('off')

plt.tight_layout()
[Link]("All_Interpolation_Comparison_1_c.png")
[Link]()
print("All three images saved in one PNG as
'All_Interpolation_Comparison.png'")

/tmp/[Link]: RuntimeWarning: invalid value

encountered in cast
bilinear_img = [Link]((R_lin, G_lin,
B_lin)).clip(0,255).astype(np.uint8)
/tmp/[Link]: RuntimeWarning: invalid value
encountered in cast
bicubic_img = [Link]((R_cub, G_cub,
B_cub)).clip(0,255).astype(np.uint8)

Image Quality Metrics (Higher Sharpness → Better Edges, Lower Noise →

Cleaner Image):
Bilinear | Sharpness: 0.027 | Noise: 14.968
Bicubic | Sharpness: 0.029 | Noise: 15.161
OpenCV Demosaic | Sharpness: 0.026 | Noise: 15.504
All three images saved in one PNG as
'All_Interpolation_Comparison.png'

OpenCV Demosaic
Aspect Bilinear Interpolation Bicubic Interpolation (RGGB)
Computati Simple, fast, uses linear More complex, uses cubic Highly optimized, very
on kernel kernel → slower fast
Complexit
y
Image Produces smoother Preserves sharper edges Balanced smoothness
Smoothne images with less noise but may cause ringing with natural color
ss artifacts rendering
Edge Slightly blurs edges Retains edges better, Best edge retention
Preservati sharper details with realistic details
on
Sharpness Lower → fewer details, Higher → retains details, Balanced → sharp
Metric fewer artifacts may amplify noise details with minimal
artifacts
Noise Handles noisy images May amplify noise if image Handles noise well due
Sensitivity better due to smoothing is noisy to internal optimization
Color No direct color No direct color correction Optimized for color
Accuracy correction fidelity in Bayer pattern
Use Case When noise reduction When sharper edges and When best visual quality
Preference and speed are preferred high-quality interpolation and real-time
are needed processing are required
''' What assumptions are you making while interpolating the missing
pixel
values and what will happen when the assumptions do not hold? '''

{"type":"string"}

Assumptions While Interpolating Missing Pixels

Smooth Color Variation: We assume that the red, green, and blue values change gradually
across neighboring pixels. This allows us to estimate missing pixels by averaging nearby known
pixels.
No Sharp Edges: The method assumes that the image does not have sudden changes in
intensity or color in a small area.

Correct CFA Pattern: We assume the Bayer pattern (e.g., RGGB) is perfectly aligned so we know
which pixels correspond to which color.

Low Noise: The observed pixel values are assumed to be accurate, with minimal sensor or
photon noise.

What Happens If Assumptions Fail

Sharp edges or textures: Bilinear or bicubic interpolation may blur edges or create color
artifacts, like “zippers” along diagonals.

Rapid color changes: Interpolated values may not match the true colors, causing unnatural
hues.

Misaligned CFA or missing pixels: Colors may be shifted or distorted.

High noise: Noise may be amplified, making the reconstructed image grainy or blotchy.

'''' e) Perform bicubic interpolation on‘[Link]’ which has been

sampled
with the CFA pattern ‘RGGB’ (provided in kodim [Link]) to reconstruct
a full
color image. '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Load [Link] and kodim_cfa.mat
# ---------------------------
raw_data = [Link]('Copy of [Link]')
cfa_data = [Link]('Copy of kodim_cfa.mat')

raw_image = raw_data['raw'].astype(np.float32) # Raw Bayer image

cfa_mask = cfa_data['cfa'] # CFA mask: 1=R,
2=G, 3=B

H, W = raw_image.shape

# ---------------------------
# Separate R, G, B channels based on CFA
# ---------------------------
R = [Link](cfa_mask == 1, raw_image, [Link])
G = [Link](cfa_mask == 2, raw_image, [Link])
B = [Link](cfa_mask == 3, raw_image, [Link])
# ---------------------------
# Coordinate grid
# ---------------------------
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

# ---------------------------
# Bicubic Interpolation Function
# ---------------------------
def bicubic_interpolate(channel):
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

# ---------------------------
# Interpolate R, G, B channels
# ---------------------------
R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)

# ---------------------------
# Reconstruct full color image
# ---------------------------
RGB_cubic = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)
RGB_cubic = RGB_cubic / [Link](RGB_cubic) # Normalize to [0,1]
RGB_cubic = np.nan_to_num(RGB_cubic) # Replace NaNs with 0

# ---------------------------
# Display and save as one PNG
# ---------------------------
fig, axes = [Link](1, 2, figsize=(12, 6))

# Raw Bayer Image

axes[0].imshow(raw_image / [Link](raw_image), cmap='gray')
axes[0].set_title('Raw Bayer Image')
axes[0].axis('off')

# Reconstructed Bicubic Image

axes[1].imshow(RGB_cubic)
axes[1].set_title('Reconstructed RGB (Bicubic Interpolation)')
axes[1].axis('off')

plt.tight_layout()
[Link]("Raw_and_Bicubic_Reconstruction.png") # Save combined
image as PNG
[Link]()

print("Image saved as 'Raw_and_Bicubic_Reconstruction_e.png'")

WARNING:[Link]:Clipping input data to the valid range for
imshow with RGB data ([0..1] for floats or [0..255] for integers). Got
range [-0.11791197982789256..1.0].

Image saved as 'Raw_and_Bicubic_Reconstruction_e.png'

''' f) Now convert the RGB image obtained in part (e) to YCrCb color
space
and median filter (choose appropriate filter size) the chrominance
channels.
Next,Convert back the image into RGB space.'''

import cv2
import numpy as np
import [Link] as plt
from [Link] import griddata
import [Link] as sio

# ---------------------------
# Load [Link] and CFA mask
# ---------------------------
raw_data = [Link]('Copy of [Link]')
cfa_data = [Link]('Copy of kodim_cfa.mat')

raw_image = raw_data['raw'].astype(np.float32)
cfa_mask = cfa_data['cfa']

H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

# ---------------------------
# Bicubic interpolation function
# ---------------------------
def bicubic_interpolate(channel):
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

# Separate R, G, B channels
R = [Link](cfa_mask == 1, raw_image, [Link])
G = [Link](cfa_mask == 2, raw_image, [Link])
B = [Link](cfa_mask == 3, raw_image, [Link])

# Bicubic interpolation
R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)

# Combine into RGB

RGB_cubic = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)
RGB_cubic = RGB_cubic / [Link](RGB_cubic)
RGB_cubic = np.nan_to_num(RGB_cubic)

# ---------------------------
# Convert to YCrCb
# ---------------------------
RGB_uint8 = (RGB_cubic * 255).astype(np.uint8)
YCrCb = [Link](RGB_uint8, cv2.COLOR_RGB2YCrCb)

# Extract channels before filtering

Y = YCrCb[:, :, 0]
Cr_before = YCrCb[:, :, 1]
Cb_before = YCrCb[:, :, 2]

# ---------------------------
# Apply median filter to Cr and Cb
# ---------------------------
filter_size = 3
Cr_after = [Link](Cr_before, filter_size)
Cb_after = [Link](Cb_before, filter_size)
# Merge back after filtering
YCrCb_filtered = [Link]([Y, Cr_after, Cb_after])
RGB_filtered = [Link](YCrCb_filtered, cv2.COLOR_YCrCb2RGB)

# ---------------------------
# Display and save as one PNG
# ---------------------------
fig, axes = [Link](2, 3, figsize=(14, 10))

# Original RGB
axes[0,0].imshow(RGB_cubic)
axes[0,0].set_title("Original Bicubic RGB")
axes[0,0].axis('off')

# Cr before
axes[0,1].imshow(Cr_before, cmap='gray')
axes[0,1].set_title("Cr Channel Before")
axes[0,1].axis('off')

# Cb before
axes[0,2].imshow(Cb_before, cmap='gray')
axes[0,2].set_title("Cb Channel Before")
axes[0,2].axis('off')

# Filtered RGB
axes[1,0].imshow(RGB_filtered)
axes[1,0].set_title("Filtered RGB")
axes[1,0].axis('off')

# Cr after
axes[1,1].imshow(Cr_after, cmap='gray')
axes[1,1].set_title("Cr Channel After")
axes[1,1].axis('off')

# Cb after
axes[1,2].imshow(Cb_after, cmap='gray')
axes[1,2].set_title("Cb Channel After")
axes[1,2].axis('off')

plt.tight_layout()
[Link]("Bicubic_YCrCb_Filtered.png") # Save entire figure as one
PNG
[Link]()

print("Figure saved as 'Bicubic_YCrCb_Filtered.png_f'")

WARNING:[Link]:Clipping input data to the valid range for

imshow with RGB data ([0..1] for floats or [0..255] for integers). Got
range [-0.11791197982789256..1.0].
Figure saved as 'Bicubic_YCrCb_Filtered.png_f'

''' g) For reference the true colour image ‘[Link]’ is given in

the data
folder. Compare the demosaiced images obtained in part (e) and (f)
with the
actual image and comment your observations. (The image ‘[Link]’
is
obtained from Kodak dataset) '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata
import cv2
from [Link] import peak_signal_noise_ratio as psnr,
structural_similarity as ssim

# ---------------------------
# Load raw Bayer data and CFA
# ---------------------------
raw_data = [Link]('Copy of [Link]')
cfa_data = [Link]('Copy of kodim_cfa.mat')

raw_image = raw_data['raw'].astype(np.float32)
cfa_mask = cfa_data['cfa']

H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

def bicubic_interpolate(channel):
"""Perform bicubic interpolation using griddata."""
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

# ---------------------------
# Separate R, G, B channels using CFA mask
# ---------------------------
R = [Link](cfa_mask == 1, raw_image, [Link])
G = [Link](cfa_mask == 2, raw_image, [Link])
B = [Link](cfa_mask == 3, raw_image, [Link])

# ---------------------------
# Bicubic Interpolation
# ---------------------------
R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)

RGB_cubic = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)

RGB_cubic = RGB_cubic / [Link](RGB_cubic)
RGB_cubic = np.nan_to_num(RGB_cubic)
RGB_cubic_uint8 = (RGB_cubic * 255).astype(np.uint8)

# ---------------------------
# Median filtering in YCrCb space
# ---------------------------
YCrCb = [Link](RGB_cubic_uint8, cv2.COLOR_RGB2YCrCb)
filter_size = 3
Y = YCrCb[:, :, 0]
Cr = [Link](YCrCb[:, :, 1], filter_size)
Cb = [Link](YCrCb[:, :, 2], filter_size)

YCrCb_filtered = [Link]([Y, Cr, Cb])

RGB_filtered = [Link](YCrCb_filtered, cv2.COLOR_YCrCb2RGB)
# ---------------------------
# Load ground-truth image (Kodak dataset)
# ---------------------------
GT = [Link]('Copy of [Link]') # use correct file from Kodak
dataset
GT = [Link](GT, cv2.COLOR_BGR2RGB)
GT = [Link](GT, (W, H)) # ensure same size for PSNR/SSIM
computation

# ---------------------------
# Compute PSNR and SSIM values
# ---------------------------
psnr_cubic = psnr(GT, RGB_cubic_uint8)
ssim_cubic = ssim(GT, RGB_cubic_uint8, channel_axis=2)

psnr_filtered = psnr(GT, RGB_filtered)

ssim_filtered = ssim(GT, RGB_filtered, channel_axis=2)

# Print comparison table

print("\nComparison with Ground Truth ([Link]):")
print(f"{'Method':<25} {'PSNR':<10} {'SSIM'}")
print("-" * 40)
print(f"{'Bicubic Interpolation':<25} {psnr_cubic:.2f}
{ssim_cubic:.3f}")
print(f"{'Chrominance Filtered':<25} {psnr_filtered:.2f}
{ssim_filtered:.3f}")

# ---------------------------
# Display and save the figure as one PNG
# ---------------------------
fig, axes = [Link](1, 3, figsize=(18, 6))

axes[0].imshow(RGB_cubic_uint8)
axes[0].set_title("Bicubic RGB")
axes[0].axis('off')

axes[1].imshow(RGB_filtered)
axes[1].set_title("Chrominance Filtered RGB")
axes[1].axis('off')

axes[2].imshow(GT)
axes[2].set_title("Ground Truth (Kodak)")
axes[2].axis('off')

plt.tight_layout()
[Link]("Bicubic_Comparison.png") # Save all 3 images in one PNG
[Link]()

print("Figure saved as 'Bicubic_Comparison_g.png'")

Comparison with Ground Truth ([Link]):
Method PSNR SSIM
----------------------------------------
Bicubic Interpolation 22.60 0.855
Chrominance Filtered 23.37 0.899

Figure saved as 'Bicubic_Comparison_g.png'

Observations

1. Bicubic Interpolation (part e): The bicubic interpolation method produces a reconstructed
full-color image that is generally smooth and visually coherent. It effectively fills in the missing
pixel values and preserves the overall structure of the image. However, edges in the image,
especially along sharp boundaries, appear slightly blurred compared to the original ground-
truth image. In areas with fine textures or high-frequency details, some color artifacts or chroma
noise are noticeable. This is because the interpolation assumes smooth color variation, which
may not hold true in textured regions.

2. Chrominance Filtered Image (part f): By applying a median filter to the chrominance channels
(Cr and Cb) in the YCrCb color space, the reconstructed image shows noticeably cleaner colors.
The filter removes small color inconsistencies and reduces chroma noise while preserving the
luminance information (Y channel), which maintains the sharpness of edges and structural
details. Compared to the bicubic-interpolated image, this approach improves both visual quality
and quantitative metrics (PSNR and SSIM). Color transitions in fine-textured regions appear
smoother, and the overall image looks closer to the true color rendition.

3. Comparison with Ground Truth ([Link]): Both reconstructed images capture the main
colors and general structure of the scene accurately. The bicubic interpolation alone is effective
but may retain minor color noise in textured regions. The chrominance-filtered image, on the
other hand, reduces these artifacts and produces a cleaner, more visually pleasing result. Edge
sharpness is better preserved, and color reproduction is more faithful to the original image.
Overall, median filtering of chrominance channels enhances the quality of the demosaiced
image without compromising the luminance detail.
# 3 White Balancing and Tone Mapping

''' a) Assume the average color of the scene to be gray and then do
white balancing.'''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Helper: Demosaicing using Bayer pattern (bicubic interpolation)
# ---------------------------
def demosaic_bayer(raw_image, bayer_mask):
"""Perform demosaicing using bicubic interpolation."""
H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

def bicubic_interpolate(channel):
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

R = [Link](bayer_mask == 1, raw_image, [Link])

G = [Link](bayer_mask == 2, raw_image, [Link])
B = [Link](bayer_mask == 3, raw_image, [Link])

R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)

RGB = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)

RGB = RGB / [Link](RGB)
RGB = np.nan_to_num(RGB)
return (RGB * 255).astype(np.uint8)

# ---------------------------
# White Balancing (Gray World Assumption)
# ---------------------------
def white_balance_gray_world(img):
"""White balance assuming average color of scene is gray."""
avgB = [Link](img[:,:,0])
avgG = [Link](img[:,:,1])
avgR = [Link](img[:,:,2])

# Mean of all channels

gray_value = (avgB + avgG + avgR) / 3
 # Scaling factors for each channel
scaleB = gray_value / avgB
scaleG = gray_value / avgG
scaleR = gray_value / avgR

# Apply scaling
img_wb = [Link]().astype(np.float32)
img_wb[:,:,0] *= scaleB
img_wb[:,:,1] *= scaleG
img_wb[:,:,2] *= scaleR

# Clip to valid range

img_wb = [Link](img_wb, 0, 255).astype(np.uint8)
return img_wb

# ---------------------------
# Load Raw Images and Bayer Masks
# ---------------------------
raw1 = [Link]("Copy of [Link]")['RawImage1']
raw2 = [Link]("Copy of [Link]")['RawImage2']
raw3 = [Link]("Copy of [Link]")['RawImage3']

bayer2 = [Link]("Copy of [Link]")['bayer2']

bayer3 = [Link]("Copy of [Link]")['bayer3']

# ---------------------------
# Step 1: Demosaic RawImage2 and RawImage3
# ---------------------------
# Convert RawImage1 to 3 channels if single-channel
RGB1 = (raw1 / [Link](raw1) * 255).astype(np.uint8)
if [Link] == 2:
RGB1 = [Link]([RGB1, RGB1, RGB1], axis=-1)

RGB2 = demosaic_bayer(raw2, bayer2)

RGB3 = demosaic_bayer(raw3, bayer3)

# ---------------------------
# Step 2: White Balancing
# ---------------------------
WB1 = white_balance_gray_world(RGB1)
WB2 = white_balance_gray_world(RGB2)
WB3 = white_balance_gray_world(RGB3)

# ---------------------------
# Step 3: Display Results
# ---------------------------
[Link](figsize=(15, 8))

# Original images
[Link](2, 3, 1); [Link](RGB1); [Link]("RawImage1
Original"); [Link]('off')
[Link](2, 3, 2); [Link](RGB2); [Link]("RawImage2
Original"); [Link]('off')
[Link](2, 3, 3); [Link](RGB3); [Link]("RawImage3
Original"); [Link]('off')

# White-balanced images
[Link](2, 3, 4); [Link](WB1); [Link]("WB RawImage1");
[Link]('off')
[Link](2, 3, 5); [Link](WB2); [Link]("WB RawImage2");
[Link]('off')
[Link](2, 3, 6); [Link](WB3); [Link]("WB RawImage3");
[Link]('off')

plt.tight_layout()
[Link]("demosaic_whitebalance_results_3_a.png", dpi=300,
bbox_inches='tight')
[Link]()

''' b) Assume that the brightest pixel to be specular highlight and

should therefore be white. Use pixel coordinate (x, y) =
i) (830, 814) for RawImage1.
ii) (1165, 280) for RawImage2.
iii) (175, 675) for RawImage3. '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Helper: Demosaicing using Bayer pattern (bicubic interpolation)
# ---------------------------
def demosaic_bayer(raw_image, bayer_mask):
H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

def bicubic_interpolate(channel):
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

R = [Link](bayer_mask == 1, raw_image, [Link])

G = [Link](bayer_mask == 2, raw_image, [Link])
B = [Link](bayer_mask == 3, raw_image, [Link])

R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)

RGB = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)

RGB = RGB / [Link](RGB)
RGB = np.nan_to_num(RGB)
return (RGB * 255).astype(np.uint8)

# ---------------------------
# White-Patch White Balancing
# ---------------------------
def white_balance_white_patch(img, bright_pixel):
r, c = bright_pixel # (row, col)
img_wb = [Link]().astype(np.float32)
scaleR = 255 / img_wb[r, c, 0]
scaleG = 255 / img_wb[r, c, 1]
scaleB = 255 / img_wb[r, c, 2]
img_wb[:,:,0] *= scaleR
img_wb[:,:,1] *= scaleG
img_wb[:,:,2] *= scaleB
return [Link](img_wb, 0, 255).astype(np.uint8)

# ---------------------------
# Load Raw Images and Bayer Masks
# ---------------------------
raw1 = [Link]("Copy of [Link]")['RawImage1']
raw2 = [Link]("Copy of [Link]")['RawImage2']
raw3 = [Link]("Copy of [Link]")['RawImage3']

bayer2 = [Link]("Copy of [Link]")['bayer2'] # grbg

bayer3 = [Link]("Copy of [Link]")['bayer3'] # rggb
# ---------------------------
# Step 1: Demosaic RawImage2 and RawImage3
# ---------------------------
RGB1 = (raw1 / [Link](raw1) * 255).astype(np.uint8)
if [Link] == 2:
RGB1 = [Link]([RGB1]*3, axis=-1)

RGB2 = demosaic_bayer(raw2, bayer2)

RGB3 = demosaic_bayer(raw3, bayer3)

# ---------------------------
# Step 2: White-Patch White Balancing
# Pixel coordinates: (x, y) -> (col, row)
WB1 = white_balance_white_patch(RGB1, (814, 830))
WB2 = white_balance_white_patch(RGB2, (280, 1165))
WB3 = white_balance_white_patch(RGB3, (675, 175))

# ---------------------------
# Step 3: Display Results
# ---------------------------
[Link](figsize=(15,8))

[Link](2,3,1)
[Link](RGB1); [Link]("RawImage1 Original"); [Link]('off')
[Link](2,3,2)
[Link](WB1); [Link]("RawImage1 White-Patch WB");
[Link]('off')

[Link](2,3,3)
[Link](RGB2); [Link]("RawImage2 Original"); [Link]('off')
[Link](2,3,4)
[Link](WB2); [Link]("RawImage2 White-Patch WB");
[Link]('off')

[Link](2,3,5)
[Link](RGB3); [Link]("RawImage3 Original"); [Link]('off')
[Link](2,3,6)
[Link](WB3); [Link]("RawImage3 White-Patch WB");
[Link]('off')

plt.tight_layout()
[Link]("white_patch_results_3_b.png", dpi=300,
bbox_inches='tight')
[Link]()
''' c) Assume that some part of the scene to be neutral. Use pixel
coordinate (x,
y) =
i) (2000, 435) for RawImage1.
ii) (445, 715) for RawImage2.
iii) (1550, 565) for RawImage3. '''

import numpy as np
import [Link] as sio
import [Link] as plt
from [Link] import griddata

# ---------------------------
# Helper: Demosaicing using Bayer pattern (bicubic interpolation)
# ---------------------------
def demosaic_bayer(raw_image, bayer_mask):
H, W = raw_image.shape
X, Y = [Link]([Link](W), [Link](H))
points = np.column_stack(([Link](), [Link]()))

def bicubic_interpolate(channel):
known_points = points[~[Link](channel).ravel()]
known_values = channel[~[Link](channel)]
return griddata(known_points, known_values, points,
method='cubic').reshape(H, W)

R = [Link](bayer_mask == 1, raw_image, [Link])

G = [Link](bayer_mask == 2, raw_image, [Link])
B = [Link](bayer_mask == 3, raw_image, [Link])

R_cubic = bicubic_interpolate(R)
G_cubic = bicubic_interpolate(G)
B_cubic = bicubic_interpolate(B)
 RGB = [Link]([R_cubic, G_cubic, B_cubic], axis=-1)
RGB = RGB / [Link](RGB)
RGB = np.nan_to_num(RGB)
return (RGB * 255).astype(np.uint8)

# ---------------------------
# Neutral-Pixel White Balancing
# ---------------------------
def white_balance_neutral_pixel(img, neutral_pixel):
"""
img: Input RGB image (H, W, 3)
neutral_pixel: Tuple (row, col) of a neutral-colored pixel
"""
r, c = neutral_pixel
img_wb = [Link]().astype(np.float32)
avg_val = [Link](img_wb[r, c, :])
scaleR = avg_val / img_wb[r, c, 0]
scaleG = avg_val / img_wb[r, c, 1]
scaleB = avg_val / img_wb[r, c, 2]
img_wb[:,:,0] *= scaleR
img_wb[:,:,1] *= scaleG
img_wb[:,:,2] *= scaleB
return [Link](img_wb, 0, 255).astype(np.uint8)

# ---------------------------
# Load Raw Images and Bayer Masks
# ---------------------------
raw1 = [Link]("Copy of [Link]")['RawImage1']
raw2 = [Link]("Copy of [Link]")['RawImage2']
raw3 = [Link]("Copy of [Link]")['RawImage3']

bayer2 = [Link]("Copy of [Link]")['bayer2'] # grbg

bayer3 = [Link]("Copy of [Link]")['bayer3'] # rggb

# ---------------------------
# Step 1: Demosaic RawImage2 and RawImage3
# ---------------------------
RGB1 = (raw1 / [Link](raw1) * 255).astype(np.uint8)
if [Link] == 2:
RGB1 = [Link]([RGB1]*3, axis=-1)

RGB2 = demosaic_bayer(raw2, bayer2)

RGB3 = demosaic_bayer(raw3, bayer3)

# ---------------------------
# Step 2: Neutral-Pixel White Balancing
# Pixel coordinates are given as (x, y) = (col, row)
WB1 = white_balance_neutral_pixel(RGB1, (435, 2000))
WB2 = white_balance_neutral_pixel(RGB2, (715, 445))
WB3 = white_balance_neutral_pixel(RGB3, (565, 1550))

# ---------------------------
# Step 3: Display Results
# ---------------------------
[Link](figsize=(15, 8))

[Link](2, 3, 1); [Link](RGB1); [Link]("RawImage1

Original"); [Link]('off')
[Link](2, 3, 2); [Link](RGB2); [Link]("RawImage2
Original"); [Link]('off')
[Link](2, 3, 3); [Link](RGB3); [Link]("RawImage3
Original"); [Link]('off')

[Link](2, 3, 4); [Link](WB1); [Link]("WB RawImage1

(Neutral Pixel)"); [Link]('off')
[Link](2, 3, 5); [Link](WB2); [Link]("WB RawImage2
(Neutral Pixel)"); [Link]('off')
[Link](2, 3, 6); [Link](WB3); [Link]("WB RawImage3
(Neutral Pixel)"); [Link]('off')

plt.tight_layout()
[Link]("neutral_pixel_whitebalance_results_3_c.png", dpi=300,
bbox_inches='tight')
[Link]()

''' d) Perform tone mapping on each of the above images using the
following two
methods:
i) Histogram equalization '''
import cv2
import numpy as np
import [Link] as plt
import [Link]

# ---------------------------
# Histogram Equalization on RGB image
# ---------------------------
def hist_equalization_rgb(img):
"""Perform histogram equalization on each RGB channel
separately."""
eq_channels = []
for i in range(3): # Loop through R, G, B channels
eq_channels.append([Link](img[:,:,i]))
return [Link](eq_channels)

# ---------------------------
# (a) Gray World White Balance
# ---------------------------
def gray_world_wb(img):
img = [Link](np.float32)
mean_b, mean_g, mean_r = [Link](img[:,:,0]), [Link](img[:,:,1]),
[Link](img[:,:,2])
mean_gray = (mean_b + mean_g + mean_r) / 3.0
scale_b, scale_g, scale_r = mean_gray / mean_b, mean_gray /
mean_g, mean_gray / mean_r

img[:,:,0] = [Link](img[:,:,0] * scale_b, 0, 255)

img[:,:,1] = [Link](img[:,:,1] * scale_g, 0, 255)
img[:,:,2] = [Link](img[:,:,2] * scale_r, 0, 255)

return [Link](np.uint8)

# ---------------------------
# (b) Specular Highlight White Balance
# ---------------------------
def specular_highlight_wb(img, coord):
img = [Link](np.float32)
x, y = coord
ref_pixel = img[y, x, :] # BGR values at given coordinate
scale = 255.0 / ref_pixel

for i in range(3):
img[:,:,i] = [Link](img[:,:,i] * scale[i], 0, 255)

return [Link](np.uint8)

# ---------------------------
# (c) Neutral Patch White Balance
# ---------------------------
def neutral_patch_wb(img, coord):
img = [Link](np.float32)
x, y = coord
ref_pixel = img[y, x, :] # BGR values at given coordinate
mean_val = [Link](ref_pixel)
scale = mean_val / ref_pixel

for i in range(3):
img[:,:,i] = [Link](img[:,:,i] * scale[i], 0, 255)

return [Link](np.uint8)

# ---------------------------
# Load RAW images from .mat files and demosaic
# ---------------------------
def load_raw_mat(path, key):
data = [Link](path)
raw = data[key]
# Normalize if 16-bit
if [Link]() > 255:
raw = (raw / [Link]() * 255).astype(np.uint8)
else:
raw = [Link](np.uint8)
# Convert Bayer to BGR (adjust pattern if needed: BG, RG, GB, GR)
rgb = [Link](raw, cv2.COLOR_BAYER_BG2BGR)
return rgb

# Load images (variable names inside .mat are RawImage1, RawImage2,

RawImage3)
Raw1 = load_raw_mat("Copy of [Link]", "RawImage1")
Raw2 = load_raw_mat("Copy of [Link]", "RawImage2")
Raw3 = load_raw_mat("Copy of [Link]", "RawImage3")

# ---------------------------
# Apply white balance methods
# ---------------------------

# (a) Gray World

WB1_a, WB2_a, WB3_a = [gray_world_wb(img) for img in [Raw1, Raw2,
Raw3]]

# (b) Specular highlight (coords given)

coords_b = [(830,814), (1165,280), (175,675)]
WB1_b = specular_highlight_wb(Raw1, coords_b[0])
WB2_b = specular_highlight_wb(Raw2, coords_b[1])
WB3_b = specular_highlight_wb(Raw3, coords_b[2])

# (c) Neutral patch (coords given)

coords_c = [(2000,435), (445,715), (1550,565)]
WB1_c = neutral_patch_wb(Raw1, coords_c[0])
WB2_c = neutral_patch_wb(Raw2, coords_c[1])
WB3_c = neutral_patch_wb(Raw3, coords_c[2])

# ---------------------------
# Apply Histogram Equalization and Display
# ---------------------------
WB_sets = [
("Gray World", [WB1_a, WB2_a, WB3_a]),
("Specular Highlight", [WB1_b, WB2_b, WB3_b]),
("Neutral Patch", [WB1_c, WB2_c, WB3_c])
]

for method_name, WB_imgs in WB_sets:

print(f"=== {method_name} White Balance + Histogram Equalization
===")
for idx, img in enumerate(WB_imgs):
hist_eq_img = hist_equalization_rgb(img)

[Link](figsize=(10,5))
[Link](1,2,1)
[Link]([Link](img, cv2.COLOR_BGR2RGB))
[Link](f"{method_name} WB Image {idx+1}")
[Link]('off')

[Link](1,2,2)
[Link]([Link](hist_eq_img, cv2.COLOR_BGR2RGB))
[Link](f"{method_name} WB + Hist Eq Image {idx+1}")
[Link]('off')

plt.tight_layout()
[Link]("tone_mapping_hist_eq_all_3_d_i.png", dpi=300,
bbox_inches='tight')
[Link]()

=== Gray World White Balance + Histogram Equalization ===

=== Specular Highlight White Balance + Histogram Equalization ===
=== Neutral Patch White Balance + Histogram Equalization ===
''' ii) Gamma correction for the gamma values of 0.5, 0.7 and 0.9 '''

import numpy as np
import [Link] as plt
import cv2

# ---------------------------
# Gamma Correction Function
# ---------------------------
def gamma_correction(img, gamma):
"""
Apply gamma correction to an RGB image.
img: Input image in uint8 format (0-255)
gamma: Gamma value
"""
img_float = img / 255.0 # Normalize to [0,1]
img_gamma = [Link](img_float, gamma)
return [Link](img_gamma * 255, 0, 255).astype(np.uint8)

# ---------------------------
# Gamma values
# ---------------------------
gamma_values = [0.5, 0.7, 0.9]

# ---------------------------
# Collect WB images from all three methods
# ---------------------------
WB_sets = {
"Gray World": [WB1_a, WB2_a, WB3_a],
"Specular Highlight": [WB1_b, WB2_b, WB3_b],
"Neutral Patch": [WB1_c, WB2_c, WB3_c]
}

# ---------------------------
# Apply gamma correction to all WB sets
# ---------------------------
for method_name, WB_imgs in WB_sets.items():
print(f"\n=== {method_name} White Balance + Gamma Correction ===")

for idx, img in enumerate(WB_imgs, start=1):

[Link](figsize=(15, 4))

for j, g in enumerate(gamma_values):
gamma_img = gamma_correction(img, g)

[Link](1, len(gamma_values), j+1)

[Link]([Link](gamma_img, cv2.COLOR_BGR2RGB))
[Link](f"Image {idx} | Gamma {g}")
[Link]('off')
 [Link](f"{method_name} WB - Image {idx}", fontsize=14)
plt.tight_layout()
[Link]("gamma_correction_all_d_ii.png", dpi=300,
bbox_inches='tight')
[Link]()

=== Gray World White Balance + Gamma Correction ===

=== Specular Highlight White Balance + Gamma Correction ===

=== Neutral Patch White Balance + Gamma Correction ===
'''' For RawImage2 use the rectangular region with coordinates: top-
left (705,924)
and bottom-right (765,984) in (row, column) order. For other raw
images choose
a region on your own and do the denoising as mentioned above.'''

import numpy as np
from [Link] import loadmat
from skimage import color, io
from tqdm import tqdm
import [Link] as plt

# ------------------------------
# Bilateral Filter Functions
# ------------------------------
def bilateral_filter_gray(img, w=5, sigma_d=3, sigma_r=0.1):
padded = [Link](img, w, mode='reflect')
filtered = np.zeros_like(img)
x, y = [Link]([Link](-w, w+1), [Link](-w, w+1))
G = [Link](-(x**2 + y**2)/(2*sigma_d**2))

for i in tqdm(range([Link][0]), desc="Bilateral Filter

(Gray)"):
for j in range([Link][1]):
local = padded[i:i+2*w+1, j:j+2*w+1]
H = [Link](-((local - img[i,j])**2)/(2*sigma_r**2))
W = G * H
filtered[i,j] = [Link](W * local) / [Link](W)
return filtered

def bilateral_filter_color(img, w=5, sigma_d=3,

sigma_r=[0.1,0.1,0.1]):
lab = color.rgb2lab(img)
padded = [Link](lab, ((w,w),(w,w),(0,0)), mode='reflect')
filtered = np.zeros_like(lab)
x, y = [Link]([Link](-w, w+1), [Link](-w, w+1))
G = [Link](-(x**2 + y**2)/(2*sigma_d**2))

for i in tqdm(range([Link][0]), desc="Bilateral Filter

(Color)"):
 for j in range([Link][1]):
local = padded[i:i+2*w+1, j:j+2*w+1, :]
filtered_pixel = []
for c in range(3):
diff = local[:,:,c] - lab[i,j,c]
H = [Link](-(diff**2)/(2*(sigma_r[c]*100)**2))
W = G * H
filtered_pixel.append([Link](W * local[:,:,c]) /
[Link](W))
filtered[i,j,:] = filtered_pixel
filtered_rgb = [Link](color.lab2rgb(filtered), 0, 1)
return filtered_rgb

def bfilter2(img, w=5, sigma=[3,0.1]):

img = [Link](np.float32)
if [Link] == 2 or ([Link]==3 and [Link][2]==1):
return bilateral_filter_gray([Link](), w=w,
sigma_d=sigma[0], sigma_r=sigma[1])
elif [Link][2] == 3:
sigma_r = sigma[1] if isinstance(sigma[1], (list, [Link]))
else [sigma[1]]*3
return bilateral_filter_color(img, w=w, sigma_d=sigma[0],
sigma_r=sigma_r)
else:
raise ValueError("Input must be grayscale or RGB image.")

# ------------------------------
# Load Raw Images
# ------------------------------
raw1 = loadmat("Copy of [Link]")
['RawImage1'].astype(np.float32)
raw2 = loadmat("Copy of [Link]")
['RawImage2'].astype(np.float32)
raw3 = loadmat("Copy of [Link]")
['RawImage3'].astype(np.float32)

raw_images = [raw1, raw2, raw3]

output_names = ["RawImage1_comparison.png",
"RawImage2_comparison.png", "RawImage3_comparison.png"]

# ------------------------------
# Define rectangular regions for noise estimation
# (row_start, col_start, row_end, col_end)
# ------------------------------
regions = [
(600, 800, 660, 860), # RawImage1 (example region)
(705, 924, 765, 984), # RawImage2 (given)
(500, 1000, 560, 1060) # RawImage3 (example region)
]
# ------------------------------
# Denoising loop (save side-by-side images)
# ------------------------------
for idx, img in enumerate(raw_images):
# Normalize to [0,1] if needed
if [Link]() > 1:
img /= 255.0

# Ensure 3rd dimension exists

if [Link] == 2:
img = img[..., [Link]]

# Extract region for noise estimation

r0, c0, r1, c1 = regions[idx]
region = img[r0:r1, c0:c1, :]
sigma_n = [Link](axis=(0,1))
sigma_r = 1.95 * sigma_n # scaling factor

# Apply bilateral filter

denoised_img = bfilter2(img, w=5, sigma=[2.5, sigma_r])

# Prepare both for saving

if denoised_img.ndim == 3 and denoised_img.shape[2] == 1:
denoised_img = denoised_img.squeeze()
img_disp = [Link]()
cmap = 'gray'
else:
img_disp = img
cmap = None

orig_uint8 = [Link](img_disp * 255, 0, 255).astype(np.uint8)

denoised_uint8 = [Link](denoised_img * 255, 0,
255).astype(np.uint8)

# Make sure both are RGB for side-by-side save

if cmap == 'gray':
orig_uint8 = [Link]([orig_uint8]*3, axis=-1)
denoised_uint8 = [Link]([denoised_uint8]*3, axis=-1)

combined = [Link]((orig_uint8, denoised_uint8))

# Save combined image as PNG

[Link](output_names[idx], combined)

# Show result
[Link](figsize=(12,5))
[Link](combined)
[Link](f"{output_names[idx]} (Original | Denoised)")
[Link]('off')
[Link]()
Bilateral Filter (Gray): 100%|██████████| 1128/1128 [01:07<00:00,
16.65it/s]

----------------------------------------------------------------------
-----
ValueError Traceback (most recent call
last)
/tmp/[Link] in <cell line: 0>()
114 denoised_uint8 = [Link]([denoised_uint8]*3, axis=-1)
115
--> 116 combined = [Link]((orig_uint8, denoised_uint8))
117
118 # Save combined image as PNG

/usr/local/lib/python3.12/dist-packages/numpy/_core/shape_base.py in
hstack(tup, dtype, casting)
356 return _nx.concatenate(arrs, 0, dtype=dtype,
casting=casting)
357 else:
--> 358 return _nx.concatenate(arrs, 1, dtype=dtype,
casting=casting)
359
360

ValueError: all the input arrays must have same number of dimensions,
but the array at index 0 has 3 dimension(s) and the array at index 1
has 2 dimension(s)

''' 1. You have to construct a Gaussian blurring kernel with a

variance σ for
denoising an image. As you know, the Gaussian function has infinite
support.
In order to implement tractability, you need to truncate the Gaussian
function
to a finite size. For a given Gaussian kernel with variance σ, at what
point will
you choose to truncate the function and why? '''

import numpy as np

def gaussian_kernel(sigma):
"""
Constructs a truncated 2D Gaussian kernel.

Args:
sigma (float): Standard deviation of the Gaussian.

Returns:
 kernel ([Link]): 2D Gaussian kernel normalized to sum=1.
"""
# Truncate at 3*sigma
radius = int(3 * sigma)
size = 2 * radius + 1 # Ensure odd size

# Create a coordinate grid

x = [Link](-radius, radius+1)
y = [Link](-radius, radius+1)
X, Y = [Link](x, y)

# Compute Gaussian function

kernel = [Link](-(X**2 + Y**2) / (2 * sigma**2))
kernel /= [Link](kernel) # Normalize to sum=1

return kernel

# Example usage:
sigma = 2.0
kernel = gaussian_kernel(sigma)
print("Gaussian Kernel:\n", kernel)
print("Kernel shape:", [Link])

Gaussian Kernel for Image Denoising

A Gaussian kernel is used in image processing for smoothing or denoising. The 2D Gaussian
function is defined as:

[ G(x, y) = \frac{1}{2 \pi \sigma^2} \exp \Big( - \frac{x^2 + y^2}{2\sigma^2} \Big) ]

• σ (sigma) is the standard deviation, controlling the spread of the Gaussian.

• Infinite support: Theoretically, the Gaussian function extends to infinity in all directions.

• Truncation: In practice, we truncate the kernel to a finite size because it is

computationally infeasible to use an infinite kernel.

Choosing the Kernel Size

A common choice is:

[ \text{Kernel size} = \lceil 6\sigma \rceil \text{ or } 2 \cdot \lceil 3\sigma \rceil + 1 ]

• Reason:
– About 99.7% of the Gaussian’s area lies within ±3σ.

– By choosing a kernel of size 6σ + 1, we include almost all significant values of

the Gaussian and ignore negligible tails.
 – Adding 1 ensures the kernel size is odd, which allows for a symmetric kernel with
a central pixel.

Key Points
1. Truncation point: ±3σ from the center captures ~99.7% of the Gaussian energy.

2. Kernel symmetry: Use an odd-sized kernel for proper alignment around the central
pixel.

3. Normalization: Always normalize so the sum of all weights equals 1, preserving the
image intensity after filtering.

4. Denoising effect: Larger σ → stronger smoothing → more noise removed, but also more
blurring.

5. Computational efficiency: Truncating reduces the number of computations compared to

using the infinite Gaussian.
''' For a given σn how would you choose the value for σr? What would
happen
if σr ≤ σn or σr≫ σn? '''

Choosing σr in Bilateral Filtering

In bilateral filtering, two parameters control the smoothing:

1. σd – spatial standard deviation (controls how much neighboring pixels contribute based
on distance)

2. σr – range standard deviation (controls how much neighboring pixels contribute based
on intensity/color difference)

Relationship Between σr and Image Noise σn

• σn: Standard deviation of the noise in the image.

• σr: Controls intensity similarity; determines how much pixels with different intensities
are averaged.

Rule of thumb:
[ \sigma_r \approx k \cdot \sigma_n ]

• Usually, ( k \in [1.5, 2] ).

• This ensures that pixels differing only due to noise are smoothed, but edges are
preserved.
Cases
1. σr ≤ σn
– The filter is too strict in intensity similarity.

– Only pixels extremely close in intensity to the center pixel are averaged.

– Effect:
• Noise is not effectively reduced.

• Edge preservation is fine but image remains noisy.

2. σr ≫ σn
– The filter considers pixels with large intensity differences as similar.

– Effect:
• Noise is reduced more aggressively.

• Edges get blurred, as pixels across edges are averaged.

• Image loses detail and sharpness.

Conclusion
• Choose σr slightly larger than σn to effectively reduce noise while preserving edges.

• Tuning σr depends on noise level and the desired trade-off between smoothing and
edge preservation.