Haralick texture extraction odd result

Good afternoon,

— Context
I am using otb within a self-made python class in order to call for various otb functions easily within a chain of operations. Among those i use the otb Haralick texture extraction which works fine. However the output seems a bit “off” (cf pictures)

Isolated pixels, but most importantly wide areas have an output value for the different haralick indices set to 0, which is the “no data value” of my file. The weirdest is that, the “patches” all have the same values across the all haralick indices. while the “isolated” pixels all behave differently.

I wonder if this is a normal behaviour of the otb haralick texture extraction.

Thank you for your help,
Jules Michard,
Regards.


—Parameters for OTB haralcik texture extraction

  • xrand and yrad 1
  • number of bins 16
  • offset default
  • haralick indices “simple”
  • min and max are set to the min/max +/- 2*std of the image

I work with OTB 7.1 and can’t update it

— Code

################## IMPORTS
import os
import time

from osgeo import gdal
from osgeo import ogr
from osgeo import osr

import numpy as np
import rasterio

##### Faire l'import d'otb à la fin sinon tout plante
import otbApplication as otb

#------------------A modifier

fname="corona_30061966"

#--------------------------- CHEMINS --------------------------
class texture:
  
  def __init__(self,filename):
    year=filename[-4:]
    self.path_corona="/home/jovyan/data_agata/"+str(year)+"/Images/"+str(filename)+"/"
    self.path_mask="/home/jovyan/data_agata/"+str(year)+"/Images/"+str(filename)+"/"  +str(filename)
    self.path_mnt='/home/jovyan/data_agata/MNT'
    self.path_indices = "/home/jovyan/data_agata/"+str(year)+"/Indices/"

    self.raster_fname=str(filename)
    self.corrected_raster=self.path_corona+str(filename)+'.tif'

#------------------ Calculs des différents indices sur l'image CORONA (composée et clippé) ------------------
  def bin_val(self,fenetre):
    if fenetre<=10:
      return int(64)
    else:
      return int(16)
  
  
#Bande: renseigner 1, 2 ,3 ou 4 pour les bandes 4,5,6, ou è respectivement
#fenetre: int

  def haralick(self,filename,bande,nom_bande,fenetre):
      
      entree= self.corrected_raster
      
      #sortie A changer pour aboutir dans dossier INDICES + ajouter taille radius
      sortie=self.path_indices + str(filename)+'_haralick_{}.tif'.format(fenetre)

      app=None
      app = otb.Registry.CreateApplication("HaralickTextureExtraction")

      #Entrée
      app.SetParameterString("in", entree)

      #Bande sur laquelle c'est calculé
      app.SetParameterInt("channel",bande)

      #Borner les indices d'Haralick (extraire max er min)
      bande_objet=rasterio.open(entree).read(bande)
      
      
      app.SetParameterInt("parameters.nbbin",self.bin_val(fenetre=fenetre))
      app.SetParameterFloat("parameters.min",float(bande_objet.mean()-2*bande_objet.std()))
      app.SetParameterFloat("parameters.max",float(bande_objet.mean()+2*bande_objet.std()))

      #Taille de la fenêtre
      app.SetParameterInt("parameters.xrad", fenetre)
      app.SetParameterInt("parameters.yrad", fenetre)
      #Offset, c'est quoi ?

      #Choix des indices  ---> Lequel choisir
      app.SetParameterString("texture","simple")

      app.SetParameterString("out", sortie)

      app.ExecuteAndWriteOutput()

      #On réaligne les étendues si décalage
      ds_ref=gdal.Open(entree)
      gt_ref=ds_ref.GetGeoTransform

      ulx, xres, xskew, uly, yskew, yres  = ds_ref.GetGeoTransform()
      lrx = ulx + (ds_ref.RasterXSize * xres)
      lry = uly + (ds_ref.RasterYSize * yres)
      
      #unset nodata ne marche pas
      commande="/opt/anaconda/envs/env_otb/bin/gdal_edit.py -a_srs epsg:32718 -a_ullr {ulx} {uly} {lrx} {lry} {sortie}".format(ulx=ulx,uly=uly,lrx=lrx,lry=lry,sortie=str(sortie))
      os.system(commande)
      ds_ref=None
    
      return
   
  ###############Execution
 
  def compute_index(self,filename):
    win_sizes=[1]
    dataset = rasterio.open(self.corrected_raster)
    

    #On a 1 bande sur CORONA
    for n_band in range(1,dataset.count+1):
      print('Pour la bande numéro {n_band}/{tot}'.format(n_band=n_band,tot=dataset.count))
      for fen in win_sizes:
  
        print("haralick pour",fen)
        self.haralick(filename,n_band,'',fen) #2 min avec hlick_rad = 11
       
     

calcul=texture(fname)
calcul.compute_index(fname)

I have managed to install a more recent version of OTB (OTB 8.0.1) in order to see if a more recent version would yiled better results.

I am using the same functions but can’t manage to run it as i encounter a sgementation fault when trying to use any “app.ExecuteWriteOutput”.

This otb version is on a docker installed from the following page : https://github.com/Terradue/docker-otb-python/tree/8.0.0

Hi @Juuuleuh

Can you try with the latest official docker image : docker pull orfeotoolbox/otb:8.1.1 ? the python bindings should work. We have fixed a segfault recently in this version

Let me know