5454from skimage .filters import sobel
5555from skimage .future import graph
5656
57+
5758#import multiprocessing as mp
5859gdal .UseExceptions ()
5960ogr .UseExceptions ()
@@ -1609,6 +1610,229 @@ def threshold(image, t):
16091610 name = 'thresholded' , colormap = 'magenta' , blending = 'additive'
16101611 )
16111612
1613+ def hough2line (inRas , outShp , edge = 'canny' , sigma = 2 ,
1614+ thresh = None , ratio = 2 , n_orient = 6 , n_scale = 5 , hArray = True , vArray = True ,
1615+ prob = False , line_length = 100 ,
1616+ line_gap = 200 , valrange = 1 , interval = 10 , band = 2 ,
1617+ min_area = None ):
1618+
1619+ # """
1620+ # Detect and write Hough lines to a line shapefile
1621+ #
1622+ # There are optionally two input arrays on the to keep line detection clean eg 2 orientations,
1623+ # such as vertical and horizontal
1624+ #
1625+ # Parameters
1626+ # ----------
1627+ #
1628+ # inRaster: string
1629+ # path to an input raster from which the geo-reffing is obtained
1630+ #
1631+ # outShp: string
1632+ # path to the output line shapefile a corresponding polygon will also be written to disk
1633+ #
1634+ # edge: string
1635+ # the edge detection method - either phase congruency or canny
1636+ # phase is default.
1637+ #
1638+ # sigma: int
1639+ # the size of stdv defining the gaussian envelope if using canny edge or phase
1640+ # a unitless value
1641+ #
1642+ # n_orient: int
1643+ # the number of orientations used if using phase congruency edge
1644+ #
1645+ # n_scale: int
1646+ # the number of scales used if using phase congruency edge
1647+ #
1648+ # vArray: bool
1649+ # whether to detect lines on approx vert axis
1650+ #
1651+ # hArray: bool
1652+ # whether to detect lines on approx horz axis
1653+ # low_t:
1654+ # the low hysteresis threshold
1655+ # the secondary low gradient threshold permitted if connected to
1656+ # a high threshold pixel
1657+ # hi_t:
1658+ # the high hysteresis threshold.
1659+ # the principal gradient threshold from which the low values are permitted
1660+ # provided they are connected to pixels of this one
1661+ #
1662+ # prob: bool
1663+ # Whether to use a probabalistic hough - default is false
1664+ #
1665+ # line_length: int
1666+ # If using prob hough the min. line length threshold
1667+ # line_gap: int
1668+ # If using prob hough the min. line gap threshold
1669+ # val_range: int
1670+ # The + - range around the orientation automatically chosen
1671+ # interval: int
1672+ # The no of intervals of the range of values tested if using auto
1673+ # """
1674+ #TODO this is FAR too long
1675+ inDataset = gdal .Open (inRas , gdal .GA_ReadOnly )
1676+
1677+ # rb = inRas.GetRasterBand(band)
1678+ rgt = inDataset .GetGeoTransform ()
1679+
1680+ pixel_res = rgt [1 ]
1681+
1682+ ref = inDataset .GetSpatialRef ()
1683+
1684+ outShapefile = outShp
1685+ outDriver = ogr .GetDriverByName ("ESRI Shapefile" )
1686+
1687+ # Remove output shapefile if it already exists
1688+ if os .path .exists (outShapefile ):
1689+ outDriver .DeleteDataSource (outShapefile )
1690+
1691+ # get the spatial ref
1692+ # ref = vlyr.GetSpatialRef()
1693+
1694+ # Create the output shapefile
1695+ outDataSource = outDriver .CreateDataSource (outShapefile )
1696+ outLayer = outDataSource .CreateLayer ("OutLyr" , geom_type = ogr .wkbMultiLineString ,
1697+ srs = ref )
1698+
1699+
1700+ # Add an ID field
1701+ idField = ogr .FieldDefn ("id" , ogr .OFTInteger )
1702+ outLayer .CreateField (idField )
1703+
1704+ empty = np .zeros ((inDataset .RasterYSize , inDataset .RasterXSize ), dtype = np .bool )
1705+
1706+ #because I am thick
1707+ #Degrees (°) Radians (rad) Radians (rad)
1708+ #0° 0 rad 0 rad
1709+ #30° π/6 rad 0.5235987756 rad
1710+ #45° π/4 rad 0.7853981634 rad
1711+ #60° π/3 rad 1.0471975512 rad
1712+ #90° π/2 rad 1.5707963268 rad
1713+ #120° 2π/3 rad 2.0943951024 rad
1714+ #135° 3π/4 rad 2.3561944902 rad
1715+ #150° 5π/6 rad 2.6179938780 rad
1716+ #180° π rad 3.1415926536 rad
1717+ #270° 3π/2 rad 4.7123889804 rad
1718+ #360° 2π rad 6.2831853072 rad
1719+
1720+
1721+ tempIm = inDataset .GetRasterBand (band ).ReadAsArray ()
1722+ bw = tempIm > 0
1723+
1724+ bwRas = inRas [:- 4 ]+ 'bw.tif'
1725+ #maskShp = inRaster[:-4]+'bwmask.shp'
1726+
1727+ #polygonize(bwRas, maskShp, outField=None, mask = True, band = 1)
1728+ props = regionprops (bw * 1 )
1729+ orient = props [0 ]['Orientation' ]
1730+
1731+ # we will need these.....
1732+ perim = mh .bwperim (bw )
1733+ # bkgrnd = invert(bw)
1734+
1735+
1736+ bw [perim == 1 ]= 0
1737+ array2raster (bw , 1 , inRas , bwRas , gdal .GDT_Byte )
1738+
1739+
1740+ # if the the binary box is pointing negatively along maj axis
1741+
1742+ if orient < 0 :
1743+ orient += np .pi
1744+
1745+ if orient < np .pi :
1746+ angleD = np .pi - orient
1747+ angleV = angleD - np .deg2rad (90 )
1748+ else :
1749+ # if the the binary box is pointing positively along maj axis
1750+ angleD = np .pi + orient
1751+ angleV = angleD + np .deg2rad (90 )
1752+
1753+
1754+
1755+ hi_t = thresh
1756+ if hi_t >= 1.0 :
1757+ low_t = np .round ((thresh / ratio ), decimals = 1 )
1758+ else :
1759+ low_t = thresh / ratio
1760+
1761+ if edge == 'phase' :
1762+ ph = _do_phasecong (tempIm , low_t , hi_t , norient = n_orient ,
1763+ nscale = n_scale , sigma = sigma )
1764+
1765+ ph [perim == 1 ]= 0
1766+
1767+ if hArray is True :
1768+ vArray = ph
1769+ if hArray is True :
1770+ hArray = ph
1771+ del ph
1772+
1773+ else :
1774+ # else it is canny
1775+ # We must have a float to get rid of zero to nonzero image
1776+ # boundary, otherwise huff will only detect the non-zero boundary
1777+ inIm = tempIm .astype (np .float32 )
1778+ inIm [inIm == 0 ]= np .nan
1779+
1780+ if hArray is True :
1781+ hArray = canny (inIm , sigma = sigma , low_threshold = low_t ,
1782+ high_threshold = hi_t )
1783+ if vArray is True :
1784+ vArray = canny (inIm , sigma = sigma , low_threshold = low_t ,
1785+ high_threshold = hi_t )
1786+ del inIm
1787+
1788+
1789+ if prob is False :
1790+ """
1791+ Standard Hough ##############################################################
1792+ """
1793+ if hasattr (vArray , 'shape' ):
1794+
1795+ empty = _std_huff (vArray , empty , outLayer , angleV , valrange , interval , rgt )#, mk=bwRas)
1796+ if hasattr (hArray , 'shape' ):
1797+ empty = _std_huff (hArray , empty , outLayer , angleD , valrange , interval , rgt )#, mk=bwRas)
1798+
1799+
1800+ else :
1801+ """
1802+ Prob Hough ##############################################################
1803+ """
1804+
1805+ if hasattr (vArray , 'shape' ):
1806+ empty = _phl_huff (vArray , empty , outLayer , angleV , valrange ,
1807+ interval , rgt , line_length , line_gap )
1808+ if hasattr (hArray , 'shape' ):
1809+ empty = _phl_huff (hArray , empty , outLayer , angleD , valrange ,
1810+ interval , rgt , line_length , line_gap )
1811+
1812+
1813+
1814+
1815+ outDataSource .SyncToDisk ()
1816+
1817+ outDataSource = None
1818+
1819+ if prob is True :
1820+ array2raster (empty , 1 , inRas , outShp [:- 3 ]+ "tif" , gdal .GDT_Int32 )
1821+ else :
1822+ inv = np .invert (empty )
1823+ inv [tempIm == 0 ]= 0
1824+ if min_area != None :
1825+ min_final = np .round (min_area / (pixel_res * pixel_res ))
1826+ if min_final <= 0 :
1827+ min_final = 4
1828+ remove_small_objects (inv , min_size = min_final , in_place = True )
1829+ #sg, _ = nd.label(inv)
1830+ segRas = outShp [:- 3 ]+ "seg.tif"
1831+ array2raster (inv , 1 , inRas , segRas , gdal .GDT_Int32 )
1832+ del tempIm , inv
1833+
1834+ polygonize (segRas , outShp [:- 4 ]+ "_poly.shp" , outField = None , mask = True , band = 1 )
1835+
16121836def colorscale (seg , prop ):
16131837
16141838 props = regionprops (np .int32 (seg ))
@@ -1896,4 +2120,4 @@ def otbMeanshift(inputImage, radius, rangeF, minSize, outShape):
18962120# os.system(cmd4)
18972121 print ('vectorisation done - process complete - phew!' )
18982122# output = subprocess.Popen([cmd], stdout=subprocess.PIPE).communicate()[0]
1899- # print(output)
2123+ # print(output)
0 commit comments