diff --git "a/D\303\251FI/DepositMethods/Camenen2013.py" "b/D\303\251FI/DepositMethods/Camenen2013.py" index de09d076e035bfc993f8f19bed0020a381a0eaee..36daf0515ca95d91ed11a227876b21e3b01971f2 100644 --- "a/D\303\251FI/DepositMethods/Camenen2013.py" +++ "b/D\303\251FI/DepositMethods/Camenen2013.py" @@ -8,7 +8,7 @@ __author__ = "Thomas DREVET" __copyright__ = "Copyright 2019, Irstea" __credits__= "Benoit CAMENEN" -__version__ = "1.0" +__version__ = "2.0" __maintainer__ = "Thomas DREVET" __email__ = "thomas.drevet@irstea.fr" __status__ = "Production" @@ -19,7 +19,7 @@ __status__ = "Production" import math import numpy as np -def Segmentation(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedimentMask): +def Segmentation(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedimentMask, HistThresholds): """ Input: + OrthoImageGS: Orthorectified image (in gray scale) - array @@ -27,15 +27,16 @@ def Segmentation(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedimentMask): + MaskSize: Size of the mask used in the method - int + Step: Size of the step used in the method - int + SedimentMask: Mask representing the sediment bed - array + + HistThresholds: vector that contains ph1, ph2, ph3 and mh1 - 4x1 vector Output: + DepositMap: Map of the sediment bed - array """ # Set the threshold values - ThreshGravel_Deposits = 155 - ThreshChannels_Deposits = 130 - ThreshChannels_Trees = 50 - LimitGravel_Deposit = 0.03 + ThreshGravel_Deposits = HistThresholds[0] + ThreshChannels_Deposits = HistThresholds[1] + ThreshChannels_Trees = HistThresholds[2] + LimitGravel_Deposit = HistThresholds[3] # Get the width and the height of the input orthorectified image width, height = OrthoImageGS.shape[:2] @@ -48,14 +49,14 @@ def Segmentation(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedimentMask): DepositMap = np.zeros((width, height)) # For all the steps in the deposit map... - for i in range (NbStepX): - for j in range (NbStepY): + for i in range (int(NbStepX)): + for j in range (int(NbStepY)): # Calculate position of the studied step Xpos = int(MaskSize/2 + (i)*Step) Ypos = int(MaskSize/2 + (j)*Step) # Extract the studied part of the orthorectified image and reshape it - WindowedImage = OrthoImageGS[Xpos-round(MaskSize/2) : Xpos+round(MaskSize/2)+1, Ypos-round(MaskSize/2) : Ypos+round(MaskSize/2)+1] + WindowedImage = OrthoImageGS[Xpos-int(round(MaskSize/2)) : Xpos+int(round(MaskSize/2))+1, Ypos-int(round(MaskSize/2)) : Ypos+int(round(MaskSize/2))+1] WindowedImage = np.reshape(WindowedImage, (-1)) # Compute the histogram and get the bins edges @@ -87,12 +88,12 @@ def Segmentation(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedimentMask): # Compute the deposit map value Val = 1*(pHist <= ThreshChannels_Trees) + \ - 40*((pHist <= ThreshChannels_Deposits)*(mHist > LimitGravel_Deposit)) + \ + 40*((pHist <= ThreshChannels_Deposits)*(pHist > ThreshChannels_Trees)*(mHist > LimitGravel_Deposit)) + \ 170*((pHist > ThreshChannels_Trees)*(pHist <= ThreshChannels_Deposits)*(mHist <= LimitGravel_Deposit)) + \ 210*((pHist > ThreshChannels_Deposits)*(pHist <= ThreshGravel_Deposits)*(mHist <= LimitGravel_Deposit)) + \ 255*((pHist > ThreshGravel_Deposits)*(mHist <= LimitGravel_Deposit)) + \ 100*((pHist > ThreshChannels_Deposits)*(mHist > LimitGravel_Deposit)) - DepositMap[round(Xpos-MaskSize/2) : round(Xpos+MaskSize/2), round(Ypos-MaskSize/2) : round(Ypos+MaskSize/2)] = Val + DepositMap[int(round(Xpos-MaskSize/2)) : int(round(Xpos+MaskSize/2)), int(round(Ypos-MaskSize/2)) : int(round(Ypos+MaskSize/2))] = Val # For all pixels from the orthorectified image... for i in range(width): diff --git "a/D\303\251FI/DepositMethods/Carbonneau2004.py" "b/D\303\251FI/DepositMethods/Carbonneau2004.py" new file mode 100644 index 0000000000000000000000000000000000000000..f16735f3f3ab2c7a6a538f26b8837adb521da7af --- /dev/null +++ "b/D\303\251FI/DepositMethods/Carbonneau2004.py" @@ -0,0 +1,103 @@ +# -*- coding: utf-8 -*- +""" +""Carbonneau2004.py: It contains the segmentation method, developped by P. Carbonneau, to detect deposits on an orthorectified image: + + + Segmentation: Method based on thresholding developped by P. Carbonneau +""" + +__author__ = "Thomas DREVET" +__copyright__ = "Copyright 2019, Irstea" +__credits__= "Patrice CARBONNEAU" +__version__ = "2.0" +__maintainer__ = "Thomas DREVET" +__email__ = "thomas.drevet@irstea.fr" +__status__ = "Production" + + +# Import the libraries +import math +import numpy as np +from skimage import color + + +def Segmentation(OrthoImage, SedimentMask, PatchSize): + """ + Input: + + OrthoImage: Orthorectified image - array + + SedimentMask: Mask representing the sediment bed - array + + PatchSize: Size of the patch used in the method - int + Output: + + DepositMap: Map of the sediment bed - array + """ + + # Convert the image in HSV + OrthoImageHSV = color.rgb2hsv(OrthoImage) + + # Extract the brightness channel + OrthoImageV = OrthoImageHSV[:, :, 2] + + # Get the width and the height of the input orthorectified image + width, height = OrthoImageV.shape[:2] + + # For all pixels from the orthorectified image... + for i in range(width): + for j in range(height): + + # Check if the pixel is inside the sediment mask... + if (SedimentMask[i, j] == 0): + # If not, Set the Bright value to null + OrthoImageV[i, j] = 0 + + # Compute the number of patches + NbPatchX = math.floor((width)/PatchSize) + NbPatchY = math.floor((height)/PatchSize) + + # Create the empty deposit map + DepositMap = np.zeros((width, height)) + + # For all the steps in the deposit map... + for PatchX in range (NbPatchX): + for PatchY in range (NbPatchY): + + # Check if the patch is inside the sediment bed + if (np.ndarray.min(SedimentMask[PatchX*PatchSize:(PatchX+1)*PatchSize, PatchY*PatchSize:(PatchY+1)*PatchSize]) == 0.0): + # if not, set the deposit map value to 0 + DepositMap[PatchX*PatchSize:(PatchX+1)*PatchSize, PatchY*PatchSize:(PatchY+1)*PatchSize] = 0.0 + # if yes... + else: + # Create an empty semi-variogram + SV = np.zeros((2*PatchSize, 2*PatchSize), dtype = float) + + # For all the lags of the semi-variogram... + for p in range (-PatchSize+1, PatchSize): + for q in range (-PatchSize+1, PatchSize): + + # Initialize the value to store the sum used in the semi-variogram + sumSV = 0 + # For all the pixels inside the patch... + for i in range (math.floor(1+(abs(p)-p)/2)+math.floor((2*PatchX)*PatchSize/2), math.floor(PatchSize-(abs(p)+p)/2+math.floor((2*PatchX)*PatchSize/2))): + for j in range (math.floor(1+(abs(q)-q)/2)+math.floor((2*PatchY)*PatchSize/2), math.floor(PatchSize-(abs(q)+q)/2)+math.floor((2*PatchY)*PatchSize/2)): + # Compute the semi-variogram sum + sumSV =+ (OrthoImageV[i+p, j+q]-OrthoImageV[i, j])**2 + # Store the semi-variogram value corresponding to the studied lag + SV[p+math.floor(PatchSize/2)-1, q+math.floor(PatchSize/2)-1] = sumSV/(2*(PatchSize-abs(p))*(PatchSize-abs(q))) + # Store in the deposit map the value corresponding to the maximum value (sill) of the semi-variogram + DepositMap[PatchX*PatchSize:(PatchX+1)*PatchSize, PatchY*PatchSize:(PatchY+1)*PatchSize] = np.ndarray.max(SV) + + # Normalize the deposit map with values between 0 and 255 + maxDM = np.max(DepositMap) + DepositMap[:, :] = np.round(DepositMap[:, :]/maxDM*255) + + # For all pixels from the orthorectified image... + for i in range(width): + for j in range(height): + + # Check if the pixel is inside the sediment mask... + if (DepositMap[i, j] == 0.0): + # Replace their value by 255 for a better display + DepositMap[i, j] = 255 + + # Convert the deposit map in int + DepositMap = DepositMap.astype(np.uint8) + # Return the deposit map + return DepositMap \ No newline at end of file diff --git "a/D\303\251FI/DepositMethods/DepositTools.py" "b/D\303\251FI/DepositMethods/DepositTools.py" new file mode 100644 index 0000000000000000000000000000000000000000..d1069ff9aca86a63424a709c3803c2d9f373b5e2 --- /dev/null +++ "b/D\303\251FI/DepositMethods/DepositTools.py" @@ -0,0 +1,118 @@ +# -*- coding: utf-8 -*- +""" +""DepositTools.py: It contains several tools that are useful to perform deposit segmentation: + + ExportSegmentAsGeoTiff: Function that permits to export a segmentation map under the Geotiff format + + GeoTiff2shpFile: Function that permits to convert a GeoTiff segmentation map into several shp files +""" + +__author__ = "Thomas DREVET" +__copyright__ = "Copyright 2019, Irstea" +__version__ = "1.5" +__maintainer__ = "Thomas DREVET" +__email__ = "thomas.drevet@irstea.fr" +__status__ = "Production" + + + +# Import the libraries +import gdal +from osgeo import osr, ogr +import numpy as np + +def ExportSegmentAsGeoTiff(SegmentMap, pathSaveSegment, MapBound, ScaleWrd2Pxl, EPSG, MaskColor): + """ + Input: + + SegmentMap: Segmentation map - array + + pathSaveSegment: path to where we have to save the map - string + + MapBound: Boundaries (xmin, xmax; ymin, ymin) of the segmentation - 2x2 array + + ScaleWrd2Pxl: Pixel size (in m) - float + + EPSG: Spatial reference - string + + MaskColor: Mask color to be removed - int + """ + # Get the width and the height of the map + width, height = SegmentMap.shape[:2] + + # Get the different classes obtained by the segmentation + Class = np.unique(SegmentMap) + + # Remove the mask color + if (MaskColor in Class): + Class = np.delete(Class, np.where(MaskColor == Class)[0][0]) + NbClass = len(Class) + + # Create a Geotiff driver in gdal, following the properties + driver = gdal.GetDriverByName("GTiff") + OutTiff = driver.Create(pathSaveSegment, height, width, (NbClass), gdal.GDT_Byte) + + # Set the Geographic transformation + OutTiff.SetGeoTransform((MapBound[0, 0], ScaleWrd2Pxl, 0, MapBound[1, 1], 0, -ScaleWrd2Pxl)) + + # Select the World Coordinate system + SpatialRef = osr.SpatialReference() + SpatialRef.ImportFromEPSG(EPSG) + wkt = SpatialRef.ExportToWkt() + OutTiff.SetProjection(wkt) + + # For all the classes (except the mask one)... + for n in range (0, NbClass): + # Initialize a empty matrix to store the binarize map of the class + BinMap = np.zeros((width, height), dtype=int) + # Compute the binarize map of the class + for i in range (width): + for j in range (height): + if (SegmentMap[i, j]==Class[n]): + BinMap[i, j] = 255 + + # Create a raster band to store the binarize map of the class + Band = OutTiff.GetRasterBand((n+1)) + # Store the binarize map of the class + Band.WriteArray(BinMap[:, :]) + # Remove the 0 values to have the transparency + Band.SetNoDataValue(0.0) + + # Clean the gdal driver + OutTiff.FlushCache() + OutTiff = None + + +def GeoTiff2shpFile(pathSegMapGeoTiff, pathSaveShp): + """ + Input: + + pathSegMapGeoTiff: path to the GeoTiff file - string + + pathSaveShp: path to where we have to save the shp file - string + """ + # Open the GeoTiff file using gdal + GeoTiffMap = gdal.Open(pathSegMapGeoTiff) + + # Extract the Spatial reference system + proj = osr.SpatialReference(wkt=GeoTiffMap.GetProjection()) + EPSG = proj.GetAttrValue('AUTHORITY', 1) + + # Extract the number of classes + NbClass = GeoTiffMap.RasterCount + + # Create a gdal driver to save the data as a shp file + Driver = ogr.GetDriverByName("ESRI Shapefile") + + # For all the classes... + for i in range (1, NbClass+1): + # Create the shp file to store the class + DataSource = Driver.CreateDataSource(pathSaveShp + "_" + str(i) + ".shp") + + # Select the World Coordinate system + SpatialRef = osr.SpatialReference() + SpatialRef.ImportFromEPSG(int(EPSG)) + + # Create the layer inside the shp file + Layer = DataSource.CreateLayer("Segmentation_" + str(i), SpatialRef) + + # Extract the data from raster band + Data = GeoTiffMap.GetRasterBand(i) + # Polygonize the data and store them in the layer + gdal.Polygonize(Data, Data, Layer, i, [], callback=None) + + # Clean the gdal driver + DataSource.FlushCache() + DataSource = None + + \ No newline at end of file diff --git "a/D\303\251FI/Orthorectification/OrthoTools.py" "b/D\303\251FI/Orthorectification/OrthoTools.py" index 5801756e43fd79a543f082b87f672d87d7cb6cf6..cd1d9b2e7ce5df171c5395772807007e13b209d1 100644 --- "a/D\303\251FI/Orthorectification/OrthoTools.py" +++ "b/D\303\251FI/Orthorectification/OrthoTools.py" @@ -5,7 +5,7 @@ + GetExtrinsicParameters: It permits to obtain the extrensic parameters from the image + FindNearest: It permits to return the value and its index of the nearest element of a matrix + ComputeOrthorectificatuion: It permits to obtain the orthorectification of the image knowing the parameters - + ExportAsGeoTiff: It permits to export the image under the geotiff format + + ExportOrthoAsGeoTiff: It permits to export the image under the geotiff format """ __author__ = "Thomas DREVET" @@ -20,8 +20,8 @@ import numpy as np import cv2 import math from ReadODS import ReadMiresCoord, ReadMiresImag -import gdal -from osgeo import osr +#import gdal +from osgeo import osr, gdal def GetExtrinsicParameters(pathMiresImageODS, pathMiresCoordODS, IntrinsicMatrix, DistortCoeffs): @@ -179,14 +179,14 @@ def ComputeOrthorectification(Image, OrthoLocation, ScaleWrd2Pxl, DTM, Intrinsic return np.transpose(OrthoImage, (1, 0, 2)) -def ExportAsGeoTiff(OrthoImage, pathSaveImage, OrthoBoundaries, ScaleWrd2Pxl, EPSG, ColorTransparency): +def ExportOrthoAsGeoTiff(OrthoImage, pathSaveImage, OrthoBoundaries, ScaleWrd2Pxl, EPSG, ColorTransparency): """ Input: + OrthoImage: Orthorectified image - array + pathSaveImage: path to where we have to save the image - string + OrthoBoundaries: Boundaries (xmin, xmax; ymin, ymin) of the orthorectification - 2x2 array + ScaleWrd2Pxl: Wanted pixel size (in m) - float - + EPSG: Spacial reference - string + + EPSG: Spacial reference - int + ColorTransparency: Color to set as transparent - 3x1 array """ @@ -203,8 +203,9 @@ def ExportAsGeoTiff(OrthoImage, pathSaveImage, OrthoBoundaries, ScaleWrd2Pxl, EP # Select the World Coordinate system SpatialRef = osr.SpatialReference() - SpatialRef.SetWellKnownGeogCS(EPSG) - OutTiff.SetProjection(SpatialRef.ExportToWkt()) + SpatialRef.ImportFromEPSG(EPSG) + wkt = SpatialRef.ExportToWkt() + OutTiff.SetProjection(wkt) # Export the R, G, B values on the raster band R = OutTiff.GetRasterBand(1) diff --git "a/D\303\251FI/Orthorectification/Orthorectification.py" "b/D\303\251FI/Orthorectification/Orthorectification.py" index 83e81cc6323eaf5b2304fad6bdbb177ec6e57fe4..ac7c7e1692564ea10b6ea76c8033451f15900904 100644 --- "a/D\303\251FI/Orthorectification/Orthorectification.py" +++ "b/D\303\251FI/Orthorectification/Orthorectification.py" @@ -17,18 +17,18 @@ import numpy as np import imageio import os import cv2 -from OrthoTools import GetExtrinsicParameters, ComputeOrthorectification, ExportAsGeoTiff +from OrthoTools import GetExtrinsicParameters, ComputeOrthorectification, ExportOrthoAsGeoTiff from ReadTXT import ReadIntParam, ReadDTM import matplotlib.pyplot as plt # Set the path to the image -pathImage = "C:\\Users\\thomas.drevet\\Documents\\Images_Cameras\\camera2\\IMAG0092.JPG" +pathImage = "C:\\Users\\thomas.drevet\\Documents\\Images_Cameras\\camera2\\IMAG0592.JPG" # Set the path to DTM file -pathDTM = "C:\\Users\\thomas.drevet\\Documents\\Seafile\\DeFI\\Données_Géoréférencement\\Flat_DTM_StMarie.txt" +pathDTM = "C:\\Users\\thomas.drevet\\Documents\\Seafile\\DeFI\\Donnees_Georeferencement\\Flat_DTM_StMarie.txt" # Set the path to the ODS file containing the mires locations -pathMiresImageODS = "C:\\Users\\thomas.drevet\\Documents\\GéoMires\\camera2\\IMAG0092.ods" +pathMiresImageODS = "C:\\Users\\thomas.drevet\\Documents\\GéoMires\\camera2\\IMAG0592.ods" # Set the path to the ODS file containing the mires world coordinates -pathMiresCoordODS = "C:\\Users\\thomas.drevet\\Documents\\Seafile\\DeFI\\Données_Géoréférencement\\coordonnees-mires-saintemarie.ods" +pathMiresCoordODS = "C:\\Users\\thomas.drevet\\Documents\\Seafile\\DeFI\\Donnees_Georeferencement\\coordonnees-mires-saintemarie.ods" # Set the pixel size (in m) ScaleWrd2Pxl = 0.05 # Set the Orthorectification boundaries @@ -104,4 +104,4 @@ if (savePNG == True): # Save the orthorectified image as a GeoTiff, if requested if (saveGeoTiff == True): - ExportAsGeoTiff(OrthoImage, "OrthoImage2_0092.tif", OrthoBoundaries, ScaleWrd2Pxl, "EPSG:3945", BkGrdColor) + ExportOrthoAsGeoTiff(OrthoImage, "OrthoImage2_0592.tif", OrthoBoundaries, ScaleWrd2Pxl, 3945, BkGrdColor) diff --git "a/D\303\251FI/TestDepositMethod.py" "b/D\303\251FI/TestDepositMethod.py" index 83c28b234e7c7309c7a8c941e3ea837cc8226d26..5c4d0523ccecd3b84c763c4e4f8c99d104ae30e1 100644 --- "a/D\303\251FI/TestDepositMethod.py" +++ "b/D\303\251FI/TestDepositMethod.py" @@ -5,7 +5,7 @@ __author__ = "Thomas DREVET" __copyright__ = "Copyright 2019, Irstea" -__version__ = "2.0" +__version__ = "2.1" __maintainer__ = "Thomas DREVET" __email__ = "thomas.drevet@irstea.fr" __status__ = "Production" @@ -16,12 +16,13 @@ __status__ = "Production" import imageio from skimage import color from DepositMethods.Camenen2013 import Segmentation as CamenenSeg -#from DepositMethods.Carbonneau2004 import Segmentation as CarbonneauSeg +from DepositMethods.Carbonneau2004 import Segmentation as CarbonneauSeg +from DepositMethods.DepositTools import ExportSegmentAsGeoTiff, GeoTiff2shpFile import matplotlib.pyplot as plt import numpy as np # Set the path to the image -pathOrthoImage = "C:\\Users\\thomas.drevet\\Documents\\Seafile\\DeFI\\Images_Orthorectifiees\\OrthoImage2_0592.png" +pathOrthoImage = "C:\\Users\\thomas.drevet\\Desktop\\OrthoImage2_0592.tif" # Set the path to the sediment mask @@ -34,28 +35,27 @@ ScaleWrd2Pxl = 0.05 MaskSize = 20 # Set the Step value -Step = 7 +Step = 10 -# Set the Patch size -PatchSize = 7 +# Set the Camenen histogram thresholds +CamenenHistThr = [155, 130, 50, 0.03] -# Set the Histogram Equalizer choice -HistStretch = True +# Set the Patch size +PatchSize = 20 # Read the orthorectified image and convert it in gray scale OrthoImage = imageio.imread(pathOrthoImage) -OrthoImageGS = color.rgb2gray(OrthoImage)*255 +OrthoImageGS = np.around(color.rgb2gray(OrthoImage)*255).astype(np.uint8) # Read the sediment mask and convert it in gray scale SedMask = imageio.imread(pathSedMask) -SedMaskGS = color.rgb2gray(SedMask)*255 +SedMaskGS = np.around(color.rgb2gray(SedMask)*255).astype(np.uint8) # Compute the deposit map with the Camenen method -DepositMapCamenen = CamenenSeg(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedMaskGS) -DepositMapCamenen.astype(np.uint8) +DepositMapCamenen = CamenenSeg(OrthoImageGS, ScaleWrd2Pxl, MaskSize, Step, SedMaskGS, CamenenHistThr) # Compute the deposit map with the Carbonneau 2004 method -DepositMapCarbonneau = CarbonneauSeg(OrthoImage, SedMaskGS, HistStretch, PatchSize) +DepositMapCarbonneau = CarbonneauSeg(OrthoImage[:, :, 0:3], SedMaskGS, PatchSize) # Display the orthorectified image and the deposit map Inp = plt.figure(0) @@ -66,10 +66,12 @@ MapCamenen = plt.figure(1) plt.imshow(DepositMapCamenen, cmap='gray') plt.title("Map Camenen") -#MapCarbonneau = plt.figure(1) -#plt.imshow(DepositMapCarbonneau , cmap='gray') -#plt.title("Map Carbonneau") +MapCarbonneau = plt.figure(2) +plt.imshow(DepositMapCarbonneau , cmap='gray') +plt.title("Map Carbonneau") # Save the deposit map -imageio.imsave("Camenen2_0592.png", DepositMapCamenen) -#imageio.imsave("Carbonneau2_0592.png", DepositMapCarbonneau) \ No newline at end of file +ExportSegmentAsGeoTiff(DepositMapCamenen, "C:\\Users\\thomas.drevet\\Documents\\Divers\\DMCamenen2_0592.tif", np.array([[1959500.00, 1959535.00], [4242750.00, 4242800.00]]), ScaleWrd2Pxl, 3945, 0) +GeoTiff2shpFile("C:\\Users\\thomas.drevet\\Documents\\Divers\\DMCamenen2_0592.tif","C:\\Users\\thomas.drevet\\Documents\\Divers\\CamenenTest") +imageio.imsave("Camenen_Drone1cm.png", DepositMapCamenen) +imageio.imsave("Carbonneau_Drone1cm_1m.png", DepositMapCarbonneau) \ No newline at end of file