Differences between revisions 3 and 19 (spanning 16 versions)
Revision 3 as of 2008-07-08 13:14:50
Size: 971
Comment:
Revision 19 as of 2008-07-11 08:17:22
Size: 7916
Comment:
Deletions are marked like this. Additions are marked like this.
Line 11: Line 11:
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. 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.
Line 27: Line 27:

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

{{{#!python
# Einlesen des Originalbildes
file='/pf/u/u242027/SAR_raw/ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
I=read_asar_imp(file)

# Darstellung des ASAR-Bildes in db (image1.png)
w=8 # Größe der "Unterboxen"
A=mean_std_box(I,w)

figure(2)
imshow(10*log10(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')
}}}

{{attachment: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.

{{{#!python
def db_hist(I,b,r):
    """Input: I ist das Intensitätsbild, b ist die Anzahl der bins, r ist die Limitierung der x-Achse"""
    A_db=nan_to_num(10*log10(I))
    h=histogram(A_db,bins=b,range=r,normed=True)
    return h[0],h[1] # h[0]: y-Achse des Histogramms, h[1]: x-Achse des Histogramms

# Definition der homogenen Bereiche
A_db=10*log10(A[:,:,0])
A1=A[295:394,929:997,0]
A2=A[198:225,168:197,0]
A3=A[863:949,382:464,0]
A4=A[846:851,623:629,0]
Aice=A[763:771,326:329,0] # vermutlich Eisberge

b=50
r=[-20,5]

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)

# Histogramm
figure(7)
plot(x1,y1,'g',x2,y2,'b',x3,y3,'y',x4,y4,'r',xice,yice,'m')
}}}

[[attachment:image2b.png]]

Die Schwellwerte haben wir aus dem Histogramm abgelesen.

{{{#!python
# Klassifikation

# Schwellwerte
s0=-17.3
s1=-12.4
s2=-7.88
s3=-1.73

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

figure(11)
imshow(A_class,cm.gist_rainbow,interpolation='nearest',origin='lower')
colorbar()
show()
}}}

[[attachment: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.

{{{#!python
from scipy import *
from pylab import *
from read_asar import *
import scipy.stats
import os.path
import scipy.ndimage as ndi


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

def smooth(I,N):
    """Box average filter"""
    kernel=ones((N,N),float32)/(N**2)
    return ndi.convolve(I, kernel)

def db_hist(I,b,r):
    A_db=nan_to_num(10*log10(I))
    h=histogram(A_db,bins=b,range=r,normed=True)
    return h[0],h[1]

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)

#*******************************************************
# Einlesen des Originalbildes

filename='/pf/u/u241110/project/asar_box_1090x1051x2.dat'
if not(os.path.exists(filename)):
    file='/pf/u/u242027/SAR_raw/ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
    I=read_asar_imp(file)
    w=8 # Größe der "Unterboxen"
    A=mean_std_box(I,w)
    A.tofile(filename)
else:
    A=reshape(fromfile(filename),(1090,1051,2))

#*********************************************************
# Darstellung des ASAR-Bildes in db (image1.png)

figure(1)
imshow(10*log10(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')
gray()
colorbar()
title('ASAR-Bild in db')

#**********************************************************
# Histogramm

#Festlegung der homogenen Flächen
A_db=db_image(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_sm=smooth(A[:,:,0],2)
A_db_sm=db_image(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]

#b=bins, r=Länge der X-Achse im Histogramm
b=50
r=[-20,5]

# Erzeugung des Histogramms
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)

# graphische Darstellung
figure(2)
plot(x1,y1,'g',x2,y2,'b',x3,y3,'y',x4,y4,'r',xice,yice,'m')
title('Histogramm (ungefiltert)')
xlabel('Intensitaet in db')
ylabel('relative Haeufigkeit')

figure(3)
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')
title('Histogramm (gefiltert)')
xlabel('Intensitaet in db')
ylabel('relative Haeufigkeit')

#*************************************************************************
# Klassifikation

# Schwellwerte
s=[-15.1,-11,-7.56,-1.98]
A_class=classification(A_db,s)

# Schwellwerte (gefiltertes Bild)
s=[-15.2,-11,-7.83,-3.09]
A_class_sm=classification(A_db_sm,s)

# graphische Darstellung
figure(4)
imshow(A_class,cm.gist_rainbow,interpolation='nearest',origin='lower')
colorbar()
title('Klassifiziertes Bild (ungefiltert)')
A_class.tofile('asar_class_1090x1051.dat')

figure(5)
imshow(A_class_sm,cm.gist_rainbow,interpolation='nearest',origin='lower')
colorbar()
title('Klassifiziertes Bild (gefiltert)')
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

{{attachment:hist_klein.png}} {{attachment:asar_class.png}}
{{attachment:hist_filt_klein.png}} {{attachment:asar_class_filt.png}}

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 ü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
  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     A_db=nan_to_num(10*log10(I))
  29     h=histogram(A_db,bins=b,range=r,normed=True)
  30     return h[0],h[1]
  31 
  32 def classification(A_db,s):
  33     """liefert klassifiziertes Bild aus db-Bild"""
  34     # s ist eine Liste mit beliebig vielen Schwellwerten
  35     A_cl=A_db.copy()
  36     A_cl[A_db<s[0]]=0
  37     for i in range(len(s)):
  38         A_cl[A_db>s[i]]=i+1
  39     return A_cl
  40 
  41 def db_image(A):
  42     """rechnet Intensität in dezibel um"""
  43     return 10*log10(A)
  44 
  45 #*******************************************************
  46 # Einlesen des Originalbildes
  47 
  48 filename='/pf/u/u241110/project/asar_box_1090x1051x2.dat'
  49 if not(os.path.exists(filename)):
  50     file='/pf/u/u242027/SAR_raw/ASA_IMP_1PNDPA20060617_043346_000000162048_00362_22460_2136.N1'
  51     I=read_asar_imp(file)
  52     w=8 # Größe der "Unterboxen"
  53     A=mean_std_box(I,w)
  54     A.tofile(filename)
  55 else:
  56     A=reshape(fromfile(filename),(1090,1051,2))
  57 
  58 #*********************************************************
  59 # Darstellung des ASAR-Bildes in db (image1.png)
  60 
  61 figure(1)
  62 imshow(10*log10(A[:,:,0]),vmin=-20,vmax=0,interpolation='nearest',origin='lower')
  63 gray()
  64 colorbar()
  65 title('ASAR-Bild in db')
  66 
  67 #**********************************************************
  68 # Histogramm
  69 
  70 #Festlegung der homogenen Flächen
  71 A_db=db_image(A[:,:,0])
  72 A1=A[365:380,976:994,0]
  73 A2=A[198:225,168:197,0]
  74 A3=A[830:860,55:80,0]
  75 A4=A[291:298,877:894,0]
  76 Aice=A[763:771,326:329,0]
  77 
  78 # homogene Flächen (gefiltertes Bild)
  79 A_sm=smooth(A[:,:,0],2)
  80 A_db_sm=db_image(A_sm)
  81 A1_sm=A_sm[365:380,976:994]
  82 A2_sm=A_sm[198:225,168:197]
  83 A3_sm=A_sm[830:860,55:80]
  84 A4_sm=A_sm[291:298,877:894]
  85 Aice_sm=A_sm[763:771,326:329]
  86 
  87 #b=bins, r=Länge der X-Achse im Histogramm
  88 b=50
  89 r=[-20,5]
  90 
  91 # Erzeugung des Histogramms
  92 y1,x1=db_hist(A1,b,r)
  93 y2,x2=db_hist(A2,b,r)
  94 y3,x3=db_hist(A3,b,r)
  95 y4,x4=db_hist(A4,b,r)
  96 yice,xice=db_hist(Aice,b,r)
  97 
  98 y1_sm,x1_sm=db_hist(A1_sm,b,r)
  99 y2_sm,x2_sm=db_hist(A2_sm,b,r)
 100 y3_sm,x3_sm=db_hist(A3_sm,b,r)
 101 y4_sm,x4_sm=db_hist(A4_sm,b,r)
 102 yice_sm,xice_sm=db_hist(Aice_sm,b,r)
 103 
 104 # graphische Darstellung
 105 figure(2)
 106 plot(x1,y1,'g',x2,y2,'b',x3,y3,'y',x4,y4,'r',xice,yice,'m')
 107 title('Histogramm (ungefiltert)')
 108 xlabel('Intensitaet in db')
 109 ylabel('relative Haeufigkeit')
 110 
 111 figure(3)
 112 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')
 113 title('Histogramm (gefiltert)')
 114 xlabel('Intensitaet in db')
 115 ylabel('relative Haeufigkeit')
 116 
 117 #*************************************************************************
 118 # Klassifikation
 119 
 120 # Schwellwerte
 121 s=[-15.1,-11,-7.56,-1.98]
 122 A_class=classification(A_db,s)
 123 
 124 # Schwellwerte (gefiltertes Bild)
 125 s=[-15.2,-11,-7.83,-3.09]
 126 A_class_sm=classification(A_db_sm,s)
 127 
 128 # graphische Darstellung
 129 figure(4)
 130 imshow(A_class,cm.gist_rainbow,interpolation='nearest',origin='lower')
 131 colorbar()
 132 title('Klassifiziertes Bild (ungefiltert)')
 133 A_class.tofile('asar_class_1090x1051.dat')
 134 
 135 figure(5)
 136 imshow(A_class_sm,cm.gist_rainbow,interpolation='nearest',origin='lower')
 137 colorbar()
 138 title('Klassifiziertes Bild (gefiltert)')
 139 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)