Differences between revisions 16 and 22 (spanning 6 versions)
Revision 16 as of 2008-07-10 12:06:26
Size: 7669
Comment:
Revision 22 as of 2008-07-11 09:58:31
Size: 8095
Comment:
Deletions are marked like this. Additions are marked like this.
Line 81: Line 81:
{{attachment:image2b.png}} [[attachment:image2b.png]]
Line 107: Line 107:
{{attachment:image3b.png}} [[attachment:image3b.png]]
Line 128: Line 128:
    """ Mittelwert über wxw Pixel"""     """ Mittelwert und Standardabweichung über wxw Pixel"""
Line 138: Line 138:
    return M     return M # liefert Mittelwert und Standardabweichung
Line 146: Line 146:
    A_db=nan_to_num(10*log10(I))     """liefert Achsen des Histogramms"""
A_db=nan_to_num(10*log10(I)) # Input nicht in db!
Line 149: Line 150:

def classification(A_db,s):
    """liefert klassifiziertes Bild aus db-Bild"""
    # s ist eine Liste mit beliebig vielen Schwellwerten
    A_cl=A_db.copy()
    A_cl[A_db<s[0]]=0
    for i in range(len(s)):
        A_cl[A_db>s[i]]=i+1
    return A_cl

def db_image(A):
    """rechnet Intensität in dezibel um"""
    return 10*log10(A)
Line 167: Line 181:
imshow(10*log10(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower') imshow(db_image(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')
Line 176: Line 190:
A_db=10*log10(A[:,:,0])
A1=A[365:380,976:994,0]
A2=
A[198:225,168:197,0]
A3=
A[830:860,55:80,0]
A4=
A[291:298,877:894,0]
Aice=
A[763:771,326:329,0]

# homogene Flächen (gefiltertes Bild)
A_db=db_image(A[:,:,0])
A_homogen=(A[365:380,976:994,0],A[198:225,168:197,0],A[830:860,55:80,0],A[291:298,877:894,0],A[763:771,326:329,0])

# gefiltertes Bild
Line 185: Line 195:
A_db_sm=10*log10(A_sm)
A1_sm=A_sm[365:380,976:994]
A2_sm=
A_sm[198:225,168:197]
A3_sm=
A_sm[830:860,55:80]
A4_sm=
A_sm[291:298,877:894]
Aice_sm=
A_sm[763:771,326:329]
A_db_sm=db_image(A_sm)
A_homogen_sm=(A_sm[365:380,976:994],A_sm[198:225,168:197],A_sm[830:860,55:80],A_sm[291:298,877:894],A_sm[763:771,326:329])
Line 197: Line 203:
y1,x1=db_hist(A1,b,r)
y2,x2=db_hist(A2,b,r)
y3,x3=db_hist(A3,b,r)
y4,x4=db_hist(A4,b,r)
yice,xice=db_hist(Aice,b,r)

y1_sm,x1_sm=db_hist(A1_sm,b,r)
y2_sm,x2_sm=db_hist(A2_sm,b,r)
y3_sm,x3_sm=db_hist(A3_sm,b,r)
y4_sm,x4_sm=db_hist(A4_sm,b,r)
yice_sm,xice_sm=db_hist(Aice_sm,b,r)
y=zeros((b,len(A_homogen)),float)
x=zeros((b,len(A_homogen)),float)
for i in range(len(A_homogen)):
    y[:,i],x[:,i]=db_hist(A_homogen[i],b,r)

y_sm=zeros((b,len(A_homogen)),float)
x_sm=zeros((b,len(A_homogen)),float)
for i in range(len(A_homogen_sm)):
    y_sm[:,i],x_sm[:,i]=db_hist(A_homogen_sm[i],b,r)
Line 211: Line 215:
plot(x1,y1,'g',x2,y2,'b',x3,y3,'y',x4,y4,'r',xice,yice,'m') plot(x[:,0],y[:,0],'g',x[:,1],y[:,1],'b',x[:,2],y[:,2],'y',x[:,3],y[:,3],'r',x[:,4],y[:,4],'m')
Line 217: Line 221:
plot(x1_sm,y1_sm,'g',x2_sm,y2_sm,'b',x3_sm,y3_sm,'y',x4_sm,y4_sm,'r',xice_sm,yice_sm,'m') plot(x_sm[:,0],y_sm[:,0],'g',x_sm[:,1],y_sm[:,1],'b',x_sm[:,2],y_sm[:,2],'y',x_sm[:,3],y_sm[:,3],'r',x_sm[:,4],y_sm[:,4],'m')
Line 226: Line 230:
s0=-15.1
s1=-11
s2=-7.56
s3=-1.98

A_class=A_db.copy()
A_class[A_db<s0]=0
A_class[A_db>s0]=1
A_class[A_db>s1]=2
A_class[A_db>s2]=3
A_class[A_db>s3]=4
s=[-15.1,-11,-7.56,-1.98]
A_class=classification(A_db,s)
Line 239: Line 234:
s0=-15.2
s1=-11
s2=-7.83
s3=-3.09

A_class_sm=A_db_sm.copy()
A_class_sm[A_db_sm<s0]=0
A_class_sm[A_db_sm>s0]=1
A_class_sm[A_db_sm>s1]=2
A_class_sm[A_db_sm>s2]=3
A_class_sm[A_db_sm>s3]=4
s=[-15.2,-11,-7.83,-3.09]
A_class_sm=classification(A_db_sm,s)
Line 265: Line 251:
Die klassifizierten Arrays sind unter /pf/u/u241110/project/asar_class_1090x1051.dat bzw. .../asar_class_filtered_1090x1051.dat

Besprechung ihrer Aufgaben ergab:

  • Einlesen als Teilaufgabe wurde von Lars geloest Gruppe arbeitet direkt mit Bildkoordinaten, welche als Arrayindex dienen

    Festlegung der Klasseneinteilung erst bei Bearbeitung der Aufgabe -> Sinnigkeitsentscheid Verwendung von Filtern und Clusterbasierte Merkmalsanalyse Ergebnis: die Zuordnung von Klassen zu Bildkoordinaten (x,y)

Dienstag

Wir haben ein Programm geschrieben, um über ein großes Bild einen Box-Filter laufen zu lassen, der Mittelwert und Standardabweichungen in einem 2-dimensionalen Array zurückgibt.

   1 def mean_std_box(I,w):
   2     """ Mittelwert über wxw Pixel"""
   3     Y,X=I.shape # Einlesen der Bilddimension
   4     M=zeros((Y/w,X/w,2),float)
   5     a=range(0,Y/w,1)
   6     b=range(0,X/w,1)
   7     for y in a:
   8         for x in b:
   9             box=I[y*w:y*w+w-1,x*w:x*w+w-1]
  10             M[y,x,0]=mean(box.flatten())
  11             M[y,x,1]=std(box.flatten())
  12     return M

Weiter haben wir uns damit beschäftigt, das Originalbild mit dem "read_asar_imp" einzulesen und in db-Darstellung auszugeben.

   1 # Einlesen des Originalbildes
   2 file='/pf/u/u242027/SAR_raw/ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
   3 I=read_asar_imp(file)
   4 
   5 # Darstellung des ASAR-Bildes in db (image1.png)
   6 w=8 # Größe der "Unterboxen"
   7 A=mean_std_box(I,w)
   8 
   9 figure(2)
  10 imshow(10*log10(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')

image1.png

Probleme hatten wir bei der Reduzierung des Speckles. Es wurde daher beschlossen diesen Schritt am Mittwoch zunächst in der ganzen Gruppe zu besprechen und dann durchzuführen.

Mittwoch

Wir haben eine Klassifizierung nur über den Mittelwert durchgeführt. Dazu haben wir das Bild in db umgerechnet und uns ein Histogramm erzeugt, bei dem jedoch keine Klassen zu unterscheiden waren. Auch eine Filterung mit dem Lee-Filter und dem Programm smooth hat am Histogramm nicht so viel verändert.

Daher haben wir aus dem gemittelten Bild homogene Bereiche ausgewählt und diese in einem Histogramm geplottet. Insgesamt haben wir fünf Klassen unterschieden, wobei die fünfte Klasse (magenta) äußerst selten vorkam, aber einen sehr hohen Rückstreukoeffizienten hatte.

   1 def db_hist(I,b,r):
   2     """Input: I ist das Intensitätsbild, b ist die Anzahl der bins, r ist die Limitierung der x-Achse"""
   3     A_db=nan_to_num(10*log10(I))
   4     h=histogram(A_db,bins=b,range=r,normed=True)
   5     return h[0],h[1] # h[0]: y-Achse des Histogramms, h[1]: x-Achse des Histogramms
   6 
   7 # Definition der homogenen Bereiche
   8 A_db=10*log10(A[:,:,0])
   9 A1=A[295:394,929:997,0]
  10 A2=A[198:225,168:197,0]
  11 A3=A[863:949,382:464,0]
  12 A4=A[846:851,623:629,0]
  13 Aice=A[763:771,326:329,0] # vermutlich Eisberge
  14 
  15 b=50
  16 r=[-20,5]
  17 
  18 y1,x1=db_hist(A1,b,r)
  19 y2,x2=db_hist(A2,b,r)
  20 y3,x3=db_hist(A3,b,r)
  21 y4,x4=db_hist(A4,b,r)
  22 yice,xice=db_hist(Aice,b,r)
  23 
  24 # Histogramm
  25 figure(7)
  26 plot(x1,y1,'g',x2,y2,'b',x3,y3,'y',x4,y4,'r',xice,yice,'m')

image2b.png

Die Schwellwerte haben wir aus dem Histogramm abgelesen.

   1 # Klassifikation
   2 
   3 # Schwellwerte
   4 s0=-17.3
   5 s1=-12.4
   6 s2=-7.88
   7 s3=-1.73
   8 
   9 A_class=A_db.copy()
  10 A_class[A_db<s0]=0
  11 A_class[A_db>s0]=1
  12 A_class[A_db>s1]=2
  13 A_class[A_db>s2]=3
  14 A_class[A_db>s3]=4
  15 
  16 figure(11)
  17 imshow(A_class,cm.gist_rainbow,interpolation='nearest',origin='lower')
  18 colorbar()
  19 show()

image3b.png

Die Klassifikation ist noch nicht perfekt.

Donnerstag

Wir haben die Klassifikation noch etwas optimiert, sodass die wichtigsten Strukturen auf dem klassifizierten Bild erkennbar sind. Wir haben außerdem noch ein Bild erstellt, das zusätzlich einen "smooth"-Filter durchlaufen hat.

Dies ist der Programmtext der vorläufig endgültigen Version.

   1 from scipy import *
   2 from pylab import *
   3 from read_asar import *
   4 import scipy.stats
   5 import os.path
   6 import scipy.ndimage as ndi
   7 
   8 
   9 def mean_std_box(I,w):
  10     """ Mittelwert und Standardabweichung über wxw Pixel"""
  11     Y,X=I.shape # Einlesen der Bilddimension
  12     M=zeros((Y/w,X/w,2),float)
  13     a=range(0,Y/w,1)
  14     b=range(0,X/w,1)
  15     for y in a:
  16         for x in b:
  17             box=I[y*w:y*w+w-1,x*w:x*w+w-1]
  18             M[y,x,0]=mean(box.flatten())
  19             M[y,x,1]=std(box.flatten())
  20     return M # liefert Mittelwert und Standardabweichung
  21 
  22 def smooth(I,N):
  23     """Box average filter"""
  24     kernel=ones((N,N),float32)/(N**2)
  25     return ndi.convolve(I, kernel)
  26 
  27 def db_hist(I,b,r):
  28     """liefert Achsen des Histogramms"""
  29     A_db=nan_to_num(10*log10(I)) # Input nicht in db!
  30     h=histogram(A_db,bins=b,range=r,normed=True)
  31     return h[0],h[1]
  32 
  33 def classification(A_db,s):
  34     """liefert klassifiziertes Bild aus db-Bild"""
  35     # s ist eine Liste mit beliebig vielen Schwellwerten
  36     A_cl=A_db.copy()
  37     A_cl[A_db<s[0]]=0
  38     for i in range(len(s)):
  39         A_cl[A_db>s[i]]=i+1
  40     return A_cl
  41 
  42 def db_image(A):
  43     """rechnet Intensität in dezibel um"""
  44     return 10*log10(A)
  45 
  46 #*******************************************************
  47 # Einlesen des Originalbildes
  48 
  49 filename='/pf/u/u241110/project/asar_box_1090x1051x2.dat'
  50 if not(os.path.exists(filename)):
  51     file='/pf/u/u242027/SAR_raw/ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
  52     I=read_asar_imp(file)
  53     w=8 # Größe der "Unterboxen"
  54     A=mean_std_box(I,w)
  55     A.tofile(filename)
  56 else:
  57     A=reshape(fromfile(filename),(1090,1051,2))
  58 
  59 #*********************************************************
  60 # Darstellung des ASAR-Bildes in db (image1.png)
  61 
  62 figure(1)
  63 imshow(db_image(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')
  64 gray()
  65 colorbar()
  66 title('ASAR-Bild in db')
  67 
  68 #**********************************************************
  69 # Histogramm
  70 
  71 #Festlegung der homogenen Flächen
  72 A_db=db_image(A[:,:,0])
  73 A_homogen=(A[365:380,976:994,0],A[198:225,168:197,0],A[830:860,55:80,0],A[291:298,877:894,0],A[763:771,326:329,0])
  74 
  75 # gefiltertes Bild
  76 A_sm=smooth(A[:,:,0],2)
  77 A_db_sm=db_image(A_sm)
  78 A_homogen_sm=(A_sm[365:380,976:994],A_sm[198:225,168:197],A_sm[830:860,55:80],A_sm[291:298,877:894],A_sm[763:771,326:329])
  79 
  80 #b=bins, r=Länge der X-Achse im Histogramm
  81 b=50
  82 r=[-20,5]
  83 
  84 # Erzeugung des Histogramms
  85 y=zeros((b,len(A_homogen)),float)
  86 x=zeros((b,len(A_homogen)),float)
  87 for i in range(len(A_homogen)):
  88     y[:,i],x[:,i]=db_hist(A_homogen[i],b,r)
  89 
  90 y_sm=zeros((b,len(A_homogen)),float)
  91 x_sm=zeros((b,len(A_homogen)),float)
  92 for i in range(len(A_homogen_sm)):
  93     y_sm[:,i],x_sm[:,i]=db_hist(A_homogen_sm[i],b,r)
  94 
  95 # graphische Darstellung
  96 figure(2)
  97 plot(x[:,0],y[:,0],'g',x[:,1],y[:,1],'b',x[:,2],y[:,2],'y',x[:,3],y[:,3],'r',x[:,4],y[:,4],'m')
  98 title('Histogramm (ungefiltert)')
  99 xlabel('Intensitaet in db')
 100 ylabel('relative Haeufigkeit')
 101 
 102 figure(3)
 103 plot(x_sm[:,0],y_sm[:,0],'g',x_sm[:,1],y_sm[:,1],'b',x_sm[:,2],y_sm[:,2],'y',x_sm[:,3],y_sm[:,3],'r',x_sm[:,4],y_sm[:,4],'m')
 104 title('Histogramm (gefiltert)')
 105 xlabel('Intensitaet in db')
 106 ylabel('relative Haeufigkeit')
 107 
 108 #*************************************************************************
 109 # Klassifikation
 110 
 111 # Schwellwerte
 112 s=[-15.1,-11,-7.56,-1.98]
 113 A_class=classification(A_db,s)
 114 
 115 # Schwellwerte (gefiltertes Bild)
 116 s=[-15.2,-11,-7.83,-3.09]
 117 A_class_sm=classification(A_db_sm,s)
 118 
 119 # graphische Darstellung
 120 figure(4)
 121 imshow(A_class,cm.gist_rainbow,interpolation='nearest',origin='lower')
 122 colorbar()
 123 title('Klassifiziertes Bild (ungefiltert)')
 124 A_class.tofile('asar_class_1090x1051.dat')
 125 
 126 figure(5)
 127 imshow(A_class_sm,cm.gist_rainbow,interpolation='nearest',origin='lower')
 128 colorbar()
 129 title('Klassifiziertes Bild (gefiltert)')
 130 A_class_sm.tofile('asar_class_filtered_1090x1051.dat')

Die klassifizierten Arrays sind unter /pf/u/u241110/project/asar_class_1090x1051.dat bzw. .../asar_class_filtered_1090x1051.dat

hist_klein.png asar_class.png hist_filt_klein.png asar_class_filt.png

LehreWiki: \AG2_ASAR_Klassifikation (last edited 2008-07-11 12:19:23 by GregorHalfmann)