Not great use of Numpy - for instance, rather than writing your coefficients out in triplicate, this needs to be a single tensor dot product. ones() is never called for in this program if you correctly use broadcasting.
Put your code in functions.
Don't hard-code system-specific file paths.
Prefer pathlib over os, for glob etc.
Don't write a (truncated) literal for pi; instead write np.pi.
Most importantly, your results are basically nonsense because the images don't only change polarisation - they change alignment. You need to apply a matching, homography and perspective de-warp to your images for this to make any sense, and OpenCV offers these to you out of the box. The demonstration code below follows https://docs.opencv.org/4.6.0/d1/de0/tutorial_py_feature_homography.html .
Suggested
from pathlib import Path
from typing import Iterable
import cv2
import matplotlib.pyplot as plt
import numpy as np
def warp_align(images: list[np.ndarray]) -> None:
sift = cv2.SIFT_create()
keys: list[tuple] = []
descriptors: list[np.ndarray] = []
for image in images:
key, descriptor = sift.detectAndCompute(image, mask=None)
keys.append(key)
descriptors.append(descriptor)
FLANN_INDEX_KDTREE = 1
flann = cv2.FlannBasedMatcher(
indexParams={'algorithm': FLANN_INDEX_KDTREE, 'trees': 5},
searchParams={'checks': 50},
)
LOWES_RATIO = 0.7
train_desc, *query_descs = descriptors
matches = [
[
m
for m, n in flann.knnMatch(query_desc, train_desc, k=2)
if m.distance < LOWES_RATIO*n.distance
]
for query_desc in query_descs
]
def keys_to_points(matched_keys: Iterable[tuple[float, float]]) -> np.ndarray:
return np.array(tuple(matched_keys), dtype=np.float32)
train_key, *query_keys = keys
for query_key, target_matches, image in zip(query_keys, matches, images[1:]):
query_points = keys_to_points(query_key[m.queryIdx].pt for m in target_matches)
train_points = keys_to_points(train_key[m.trainIdx].pt for m in target_matches)
M, mask = cv2.findHomography(query_points, train_points, method=cv2.RANSAC, ransacReprojThreshold=5)
print(M, '\n')
cv2.warpPerspective(src=image, dst=image, M=M, dsize=image.shape[1::-1])
def calculate_polarimetry(images: list[np.ndarray]) -> np.ndarray:
all_images = np.stack(images)
coefficients = np.array((0.30, 0.59, 0.11))
i00, i45, i90 = np.tensordot(all_images, coefficients, axes=[3, 0])
stokes_i = i00 + i90
stokes_q = i00 - i90
stokes_u = 2*i45 - stokes_i
polint = np.sqrt(stokes_q*stokes_q + stokes_u*stokes_u)
poldolp = polint/(stokes_i + 0.001)
polaop = 0.5 * np.arctan(stokes_u, stokes_q)
h = polaop/np.pi + 0.5
s = poldolp
v = polint
hsv_polar = 255*np.clip(cv2.merge((h, s, v)), 0, 1)
return cv2.cvtColor(hsv_polar.astype('uint8'), cv2.COLOR_HSV2RGB)
def main() -> None:
paths = sorted(Path.cwd().glob("pol*.png"))
images = [
cv2.cvtColor(cv2.imread(str(path)), cv2.COLOR_BGR2RGB)
for path in paths
]
warp_align(images)
rgb_img = calculate_polarimetry(images)
fig, ((tl, tr), (bl, br)) = plt.subplots(nrows=2, ncols=2)
tl.imshow(images[0])
tr.imshow(images[1])
bl.imshow(images[2])
br.imshow(rgb_img)
plt.show()
# cv2.imwrite("hsvpolarized.jpg", rgbimg)
if __name__ == '__main__':
main()
Prior to homographic correction
After homographic correction
