Initial Commit

requirements.txt

@@ -0,0 +1,7 @@
+numpy
+pillow
+cython
+matplotlib
+opencv-python
+scipy
+shapely

rooftop.py

@@ -0,0 +1,320 @@
+from __future__ import print_function
+# import Image
+from PIL import Image
+import cv2
+import matplotlib.pyplot as plt
+import numpy as np
+import math 
+import glob
+from shapely.geometry import Polygon
+
+zoom = 20
+tileSize = 256
+initialResolution = 2 * math.pi * 6378137 / tileSize
+originShift = 2 * math.pi * 6378137 / 2.0
+earthc = 6378137 * 2 * math.pi
+factor = math.pow(2, zoom)
+map_width = 256 * (2 ** zoom)
+
+
+def grays(im):
+ return cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
+
+
+def white_image(im):
+ return cv2.bitwise_not(np.zeros(im.shape, np.uint8))
+
+
+def pixels_per_mm(lat, length):
+ return length / math.cos(lat * math.pi / 180) * earthc * 1000 / map_width
+
+
+def sharp(gray):
+ blur = cv2.bilateralFilter(gray, 5, sigmaColor=7, sigmaSpace=5)
+ kernel_sharp = np.array((
+ [-2, -2, -2],
+ [-2, 17, -2],
+ [-2, -2, -2]), dtype='int')
+ return cv2.filter2D(blur, -1, kernel_sharp)
+
+
+def solar_panel_params():
+ panel_lens = input("Number of panels together: ")
+ # panel_wids = input("Number of panels in breadth: ")
+ panel_wids = 1
+ length_s = input("Enter length of panel in mm: ")
+ width = input("Enter width of panel in mm: ")
+ angle = input("Rotation Angle for Solar Panels: ")
+ return panel_lens, panel_wids, length_s, width, angle
+
+
+def contours_canny(cnts):
+ cv2.drawContours(canny_contours, cnts, -1, 255, 1)
+
+ # Removing the contours detected inside the roof
+ for cnt in cnts:
+ counters = 0
+ cnt = np.array(cnt)
+ cnt = np.reshape(cnt, (cnt.shape[0], cnt.shape[2]))
+ pts = []
+
+ if cv2.contourArea(cnt) > 10:
+ for i in cnt:
+ x, y = i
+ if edged[y, x] == 255:
+ counters += 1
+ pts.append((x, y))
+
+ if counters > 10:
+ pts = np.array(pts)
+ pts = pts.reshape(-1, 1, 2)
+ cv2.polylines(canny_polygons, [pts], True, 0)
+
+
+def contours_img(cnts):
+ cv2.drawContours(image_contours, cnts, -1, 255, 1)
+
+ # Removing the contours detected inside the roof
+ for cnt in cnts:
+ counter = 0
+ cnt = np.array(cnt)
+ cnt = np.reshape(cnt, (cnt.shape[0], cnt.shape[2]))
+ pts = []
+ if cv2.contourArea(cnt) > 5:
+ for i in cnt:
+ x, y = i
+ if edged[y, x] == 255:
+ counter += 1
+ pts.append((x, y))
+ if counter > 10:
+ pts = np.array(pts)
+ pts = pts.reshape(-1, 1, 2)
+ cv2.polylines(image_polygons, [pts], True, 0)
+
+
+def rotation(center_x, center_y, points, ang):
+ angle = ang * math.pi / 180
+ rotated_points = []
+ for p in points:
+ x, y = p
+ x, y = x - center_x, y - center_y
+ x, y = (x * math.cos(angle) - y * math.sin(angle), x * math.sin(angle) + y * math.cos(angle))
+ x, y = x + center_x, y + center_y
+ rotated_points.append((x, y))
+ return rotated_points
+
+
+def createLineIterator(P1, P2, img):
+ imageH = img.shape[0]
+ imageW = img.shape[1]
+ P1X = P1[0]
+ P1Y = P1[1]
+ P2X = P2[0]
+ P2Y = P2[1]
+
+ # difference and absolute difference between points
+ # used to calculate slope and relative location between points
+ dX = P2X - P1X
+ dY = P2Y - P1Y
+ dXa = np.abs(dX)
+ dYa = np.abs(dY)
+
+ # predefine numpy array for output based on distance between points
+ itbuffer = np.empty(shape=(np.maximum(dYa, dXa), 3), dtype=np.float32)
+ itbuffer.fill(np.nan)
+
+ # Obtain coordinates along the line using a form of Bresenham's algorithm
+ negY = P1Y > P2Y
+ negX = P1X > P2X
+ if P1X == P2X: # vertical line segment
+ itbuffer[:, 0] = P1X
+ if negY:
+ itbuffer[:, 1] = np.arange(P1Y - 1, P1Y - dYa - 1, -1)
+ else:
+ itbuffer[:, 1] = np.arange(P1Y + 1, P1Y + dYa + 1)
+ elif P1Y == P2Y: # horizontal line segment
+ itbuffer[:, 1] = P1Y
+ if negX:
+ itbuffer[:, 0] = np.arange(P1X - 1, P1X - dXa - 1, -1)
+ else:
+ itbuffer[:, 0] = np.arange(P1X + 1, P1X + dXa + 1)
+ else: # diagonal line segment
+ steepSlope = dYa > dXa
+ if steepSlope:
+ slope = dX.astype(float) / dY.astype(float)
+ if negY:
+ itbuffer[:, 1] = np.arange(P1Y - 1, P1Y - dYa - 1, -1)
+ else:
+ itbuffer[:, 1] = np.arange(P1Y + 1, P1Y + dYa + 1)
+ itbuffer[:, 0] = (slope * (itbuffer[:, 1] - P1Y)).astype(int) + P1X
+ else:
+ slope = dY.astype(float) / dX.astype(float)
+ if negX:
+ itbuffer[:, 0] = np.arange(P1X - 1, P1X - dXa - 1, -1)
+ else:
+ itbuffer[:, 0] = np.arange(P1X + 1, P1X + dXa + 1)
+ itbuffer[:, 1] = (slope * (itbuffer[:, 0] - P1X)).astype(int) + P1Y
+
+ # Remove points outside of image
+ colX = itbuffer[:, 0]
+ colY = itbuffer[:, 1]
+ itbuffer = itbuffer[(colX >= 0) & (colY >= 0) & (colX < imageW) & (colY < imageH)]
+
+ # Get intensities from img ndarray
+ itbuffer[:, 2] = img[itbuffer[:, 1].astype(np.uint), itbuffer[:, 0].astype(np.uint)]
+
+ return itbuffer
+
+
+def panel_rotation(panels_series, solar_roof_area):
+
+ high_reso = cv2.pyrUp(solar_roof_area)
+ rows, cols = high_reso.shape
+ high_reso_new = cv2.pyrUp(new_image)
+
+ for _ in range(panels_series - 2):
+ for col in range(0, cols, l + 1):
+ for row in range(0, rows, w + 1):
+
+ # Rectangular Region of interest for solar panel area
+ solar_patch = high_reso[row:row + (w + 1) * pw + 1, col:col + ((l * pl) + 3)]
+ r, c = solar_patch.shape
+
+ # Rotation of rectangular patch according to the angle provided
+ patch_rotate = np.array([[col, row], [c + col, row], [c + col, r + row], [col, r + row]], np.int32)
+ rotated_patch_points = rotation((col + c) / 2, row + r / 2, patch_rotate, solar_angle)
+ rotated_patch_points = np.array(rotated_patch_points, np.int32)
+
+ # Check for if rotated points go outside of the image
+ if (rotated_patch_points > 0).all():
+ solar_polygon = Polygon(rotated_patch_points)
+ polygon_points = np.array(solar_polygon.exterior.coords, np.int32)
+
+ # Appending points of the image inside the solar area to check the intensity
+ patch_intensity_check = []
+
+ # Point polygon test for each rotated solar patch area
+ for j in range(rows):
+ for k in range(cols):
+ if cv2.pointPolygonTest(polygon_points, (k, j), False) == 1:
+ patch_intensity_check.append(high_reso[j, k])
+
+ # Check for the region available for Solar Panels
+ if np.mean(patch_intensity_check) == 255:
+
+ # Moving along the length of line to segment solar panels in the patch
+ solar_line_1 = createLineIterator(rotated_patch_points[0], rotated_patch_points[1], high_reso)
+ solar_line_1 = solar_line_1.astype(int)
+ solar_line_2 = createLineIterator(rotated_patch_points[3], rotated_patch_points[2], high_reso)
+ solar_line_2 = solar_line_2.astype(int)
+ line1_points = []
+ line2_points = []
+ if len(solar_line_2) > 10 and len(solar_line_1) > 10:
+
+ # Remove small unwanted patches
+ cv2.fillPoly(high_reso, [rotated_patch_points], 0)
+ cv2.fillPoly(high_reso_new, [rotated_patch_points], 0)
+ cv2.polylines(high_reso_orig, [rotated_patch_points], 1, 0, 2)
+ cv2.polylines(high_reso_new, [rotated_patch_points], 1, 0, 2)
+
+ cv2.fillPoly(high_reso_orig, [rotated_patch_points], (0, 0, 255))
+ cv2.fillPoly(high_reso_new, [rotated_patch_points], (0, 0, 255))
+
+ for i in range(5, len(solar_line_1), 5):
+ line1_points.append(solar_line_1[i])
+ for i in range(5, len(solar_line_2), 5):
+ line2_points.append(solar_line_2[i])
+
+ # Segmenting Solar Panels in the Solar Patch
+ for points1, points2 in zip(line1_points, line2_points):
+ x1, y1, _ = points1
+ x2, y2, _ = points2
+ cv2.line(high_reso_orig, (x1, y1), (x2, y2), (0, 0, 0), 1)
+ cv2.line(high_reso_new, (x1, y1), (x2, y2), (0, 0, 0), 1)
+
+ # Number of Solar Panels in series (3/4/5)
+ panels_series = panels_series - 1
+ result = Image.fromarray(high_reso_orig)
+ resut_2 = Image.fromarray(high_reso_new)
+ result.save('output' + fname )
+ resut_2.save('panels' + fname)
+ plt.figure()
+ plt.axis('off')
+ plt.imshow(high_reso_orig)
+ plt.figure()
+ plt.axis('off')
+ plt.imshow(high_reso_new)
+ plt.show()
+
+
+if __name__ == "__main__":
+ images = glob.glob('1.jpg')
+ # latitude = ??
+ # pl, pw, l, w, solar_angle = solar_panel_params()
+ # length, width = pixels_per_mm(latitude)
+
+ for fname in images:
+ # pl = No of panels together as length commonside, pw = Same as for pw here w = width
+ # l = Length of panel in mm, w = Width of panel in mm
+ # solar_angle = Angle for rotation
+ pl, pw, l, w, solar_angle = 4, 1, 8, 5, 30
+ image = cv2.imread(fname)
+ img = cv2.pyrDown(image)
+ print('image shape : ',img.shape)
+ n_white_pix = np.sum(img==255)
+ print('num of white pixels : ',n_white_pix)
+ # Upscaling of Image
+ high_reso_orig = cv2.pyrUp(image)
+
+ # White blank image for contours of Canny Edge Image
+ canny_contours = white_image(image)
+ # White blank image for contours of original image
+ image_contours = white_image(image)
+
+ # White blank images removing rooftop's obstruction
+ image_polygons = grays(canny_contours)
+ canny_polygons = grays(canny_contours)
+
+ # Gray Image
+ grayscale = grays(image)
+ plt.figure()
+ plt.title('grayscale')
+ plt.imshow(image, cmap='gray')
+ # Edge Sharpened Image
+ sharp_image = sharp(grayscale)
+ plt.figure()
+ plt.title('sharp_image')
+ plt.imshow(sharp_image, cmap='gray')
+
+ # Canny Edge
+ edged = cv2.Canny(sharp_image, 180, 240)
+ plt.figure()
+ plt.title('edge_image')
+ plt.imshow(sharp_image, cmap='gray') 
+ # Otsu Threshold (Automatic Threshold)
+ thresh = cv2.threshold(sharp_image, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]
+ plt.figure()
+ plt.title('Threshold_image')
+ plt.imshow(sharp_image, cmap='gray')
+ # Contours in Original Image
+ contours_img(cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)[-2])
+ # Contours in Canny Edge Image
+ contours_canny(cv2.findContours(edged, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)[-2])
+
+ # Optimum place for placing Solar Panels
+ solar_roof = cv2.bitwise_and(image_polygons, canny_polygons)
+ plt.figure()
+ plt.title('solar_roof_area')
+ plt.imshow(solar_roof, cmap='gray')
+ n_white_pix = np.sum(solar_roof==255)
+ print('solar white pix : ',n_white_pix)
+ print('size of solar roof : ',solar_roof.shape)
+ new_image = white_image(image)
+ plt.figure()
+ plt.title('new_image')
+ plt.imshow(new_image, cmap='gray')
+ print('new image shape',new_image.shape)
+ # Rotation of Solar Panels
+ panel_rotation(pl, solar_roof)
+
+ plt.show()