Preparazione dati per segmentazione rete neurale

 Ho trovato in rete un dataset di foto da drone con gia' classificate le fratture. Si tratta di dati con licenza permissiva 

 

 

I dati derivano da questo articolo 

https://www.researchgate.net/publication/357994422_A_new_subsampling_methodology_to_optimize_the_characterization_of_two-dimensional_bedrock_fracture_networks

@article{article,
author = {Ovaskainen, Nikolas and Nordbäck, Nicklas and Skyttä, Pietari and Engström, Jon},
year = {2022},
month = {01},
pages = {104528},
title = {A new subsampling methodology to optimize the characterization of two-dimensional bedrock fracture networks},
volume = {155},
journal = {Journal of Structural Geology},
doi = {10.1016/j.jsg.2022.104528}
}

e sono disponibili come ortofoto RGB da drone a questo link https://zenodo.org/records/4719627  (16 Gb) mentre le fratture sono tematizzate in un diversi file gpkg a questo link https://www.kaggle.com/datasets/nasurz/getaberget-fracture-trace-dataset/data (4Mb)

il problema e' che in questo formato le immagini non sono adatte al processamento per la rete neurale. L'idea e' quindi di dividere le ortofoto in tasselli da 15 metri e creare una immagine maschera con le sole fratturazioni. Tale compito e' stato eseguito con uno script Python in QGis

Nello script si devono impostare le coordinate geografiche della finestra  ed il folder dove saranno salvate le varie tiles. In un primo passaggio si disattiva dal progetto QGis il tema lineare delle fratture per avere solo i tasselli delle ortofoto e dopo si disattiva il tematismo raster per creare le maschere

 ------------------------------------------------------------------------------------------

from qgis.core import (
    QgsProject, QgsApplication, QgsLayout, QgsPrintLayout, QgsLayoutItemMap,
    QgsUnitTypes, QgsLayoutSize, QgsLayoutPoint, QgsLayoutExporter, QgsRectangle
)
import os



# Path to your QGIS project
project_path = "C:\\Users\\l.innocenti\\Desktop\\geta.qgz"

# Output folder
output_dir = "C:\\Users\\l.innocenti\\Desktop\\geta5\\mask\\"
os.makedirs(output_dir, exist_ok=True)

# Load project
project = QgsProject.instance()
project.read(project_path)


from qgis.core import QgsRectangle

# Compute project extent from all layers
extent = None
for layer in project.mapLayers().values():
    if not layer.isValid():
        continue
    layer_extent = layer.extent()
    if extent is None:
        extent = QgsRectangle(layer_extent)
    else:
        extent.combineExtentWith(layer_extent)

if extent is None:
    raise ValueError("No valid layers found in the project.")

xmin, ymin = 108914, 6720011
xmax, ymax = 108969, 6720066

# Get the full extent in project CRS
#extent = project.boundingBox()
#xmin, ymin, xmax, ymax = extent.xMinimum(), extent.yMinimum(), extent.xMaximum(), extent.yMaximum()

# Tile size in meters
tile_size = 15

# Iterate over the grid
x = xmin
tile_id = 0
while x < xmax:
    y = ymin
    while y < ymax:
        # Tile extent
        tile_extent = QgsRectangle(x, y, x + tile_size, y + tile_size)

        # Create a print layout
        layout = QgsPrintLayout(project)
        layout.initializeDefaults()
        layout.setName(f"tile_{tile_id}")
        pc = layout.pageCollection().pages()[0]

        # Map item
        map_item = QgsLayoutItemMap(layout)
        map_item.setRect(20, 20, 200, 200)
        map_item.setExtent(tile_extent)
        map_item.setBackgroundColor(QColor(255, 255, 255, 0))
        map_item.attemptResize(QgsLayoutSize(200, 200, QgsUnitTypes.LayoutMillimeters))
        visible_layers = [lyr for lyr in iface.mapCanvas().layers()]
        map_item.setLayers(visible_layers)

        layout.addLayoutItem(map_item)
        layout.addLayoutItem(map_item)
        map_item.attemptMove(QgsLayoutPoint(0, 0, QgsUnitTypes.LayoutMillimeters))

        # Export PNG
        exporter = QgsLayoutExporter(layout)
        out_path = os.path.join(output_dir, f"tile_{tile_id}.png")
        exporter.exportToImage(out_path, QgsLayoutExporter.ImageExportSettings())

        print(f"Saved {out_path}")
        tile_id += 1

        y += tile_size
    x += tile_size
  ------------------------------------------------------------------------------------------

 alla fine si hanno coppie di questo tipo (le immagini sono state rifilate con imagemagick per togliere lo spazio bianco inserito da QGis)

 

 

 

 per "aggiustare" il background con colore nero e il target come grigio come 127 

import cv2
import numpy as np
import sys

# --- Check for input filename ---
if len(sys.argv) < 2:
    print("Usage: python script.py <image_filename>")
    sys.exit(1)

filename = sys.argv[1]

# Read image
img = cv2.imread(filename)

if img is None:
    raise FileNotFoundError(f"Cannot read file: {filename}")

# Define white color in BGR
white = np.array([255, 255, 255], dtype=np.uint8)

# Create mask for white pixels
mask_white = np.all(img == white, axis=-1)

# Create output image
output = np.full_like(img, (128, 128, 128))  # Set all to gray
output[mask_white] = (0, 0, 0)  # Set white pixels to black

# Save with same filename (overwrite original)
cv2.imwrite(filename, output)

print(f"Image processed and saved: {filename}")