# Last edited on 2012-01-04 22:06:10 by stolfilocal

Obtaining decent satellite images of Belo Monte region

BASICS

  tdir="${HOME}/projects/img-belo-monte"

FETCHING IMAGES

  Looked in INPE's catalog for images of the region:
  
    http://www.dgi.inpe.br/CDSR/
 
  Fetched from INPE the following images of the region
  (see "wget-images.sh")
  
    CBERS_2B_HRC_20090807_164_B_103_3
    CBERS_2B_HRC_20090807_164_B_104_1
    
    LANDSAT_5_TM_20110704_225_062
    LANDSAT_5_TM_20110704_225_063

    LANDSAT_5_TM_20110727_226_062
    LANDSAT_5_TM_20110727_226_063
    
    LANDSAT_5_TM_20110913_226_062
  
  The CBERS images are large (13082 x 11536) and higher-resolution
  but do not quite cover the desired area. Deleted them.
  
  The LANDSAT images are smaller (about 7900x6700) 
  but cover neatly the desired area. 
  
  The LANDSAT_5_TM_20110913_226_062 image covers the same area
  as LANDSAT_5_TM_20110727_226_062 but about 2 months later.
  It has more clouds.  Deleted it.
  
  Created short links to the remaining four images (see "make-links.sh"):

    L52a LANDSAT_5_TM_20110704_225_062 7954 x 6739
    L62a LANDSAT_5_TM_20110727_226_062 7858 x 6723
    L53a LANDSAT_5_TM_20110704_225_063 7964 x 6704
    L63a LANDSAT_5_TM_20110727_226_063 7989 x 6711
    
  The BAND6 images are even smaller (about 1/4 linear size).

CONVERTING AND SCALING IMAGES

  Converted all images from TIFF format to PNG format using
  ImageMagick/convert. Also produced smaller versions (50%, 25%, 5%).
  See "reduce-images.sh".
  
TEXTURE AND VALUE ANALYSIS

  Two files of representative points in L225a were gathered, for forest
  and farmland. Each point is the center of a 32 x 32 pixel
  window of the same soil type.
  
    L52a/samples-trees.txt
    L52a/samples-farms.txt

HISTOGRAMS

  Computing histograms of the raw landsat images, top and bottom halves:

    cd ${tdir}

    bands=( 1 2 3 4 5 6 7 )
    bcolors=( '#ff0000' '#bb4400' '#447700' '#009900' '#007777' '#0033ff' '#8800ff' '#ff0099' )
    for band in ${bands[@]} ; do
      i=0
      if [[ "/${band}" == "/6" ]]; then iwd=490; hht=0200; else iwd=1900; hht=0800; fi
      plotcmd="plot (1.00) with lines lc rgb \"#555555\"";
      for img in L52a L62a L53a L63a ; do 
        for top in '0000' "${hht}"; do
          hfile="${img}/t${band}+${top}.hist"
          pngtopnm ${img}/L2_BAND${band}-r025.png \
            | pnmdepth 255 \
            | pnmcut 0 ${top} ${iwd} ${hht} \
            | pgmhist \
            | egrep -e '^[ ]*[0-9]' \
            | tr -d '%' \
            > ${hfile}
            plotcmd="${plotcmd}, \"${hfile}\" using 1:3 title \"${img} ${top}\" with histeps lt 1 lc rgb \"${bcolors[${i}]}\""
            i=$(( ${i} + 1 ))
        done
      done
      gnuplot \
        -e "set size 1,0.75" \
        -e "set title \"band ${band}\"" \
        -e "set logscale y" \
        -e "set xtics 20; set mxtics 2; set grid xtics mxtics" \
        -e "${plotcmd}" \
        -e "pause 300" &
    done    

  Checking the safe vmin,vmax:

    cd ${tdir}
    
    bands=( 1 2 3 4 5 6 7 )
    
    declare -a vmin
    vmin[1]=40
    vmin[2]=10
    vmin[3]=0
    vmin[4]=0
    vmin[5]=0
    vmin[6]=115
    vmin[7]=0
    
    declare -a vmax
    vmax[1]=70
    vmax[2]=42
    vmax[3]=42
    vmax[4]=120
    vmax[5]=110
    vmax[6]=144
    vmax[7]=49

    for img in L52a L62a L53a L63a ; do 
      for band in ${bands[@]} ; do
        tmin="`echo ${vmin[${band}]}/255 | bc -lq`"
        tmax="`echo ${vmax[${band}]}/255 | bc -lq`"
        echo "image = ${img} band = ${band}"
        convert ${img}/L2_BAND${band}-r025.png PGM:- \
          | pnmdepth 255 \
          | pamthreshold -simple -threshold="${tmin}" \
          | pamtopnm \
          | pnmdepth 255 \
          > ${img}/t${band}-th-min.pgm
        convert ${img}/L2_BAND${band}-r025.png PGM:- \
          | pnmdepth 255 \
          | pamthreshold -simple -threshold="${tmax}" \
          | pamtopnm \
          | pnmdepth 255 \
          > ${img}/t${band}-th-max.pgm
        pnmxarith \
          -add -scale 0.5 \
          ${img}/t${band}-th-min.pgm \
          ${img}/t${band}-th-max.pgm \
          > ${img}/t${band}-th.pgm
        pgmhist ${img}/t${band}-th.pgm
        convert ${img}/L2_BAND${band}-r025.png PGM:- \
          | pnmdepth 255 \
          | pnmnorm \
              -bvalue=${vmin[${band}]} \
              -wvalue=${vmax[${band}]} \
          > ${img}/t${band}-tn.pgm
        display -resize '25%' ${img}/t${band}-th.pgm ${img}/t${band}-tn.pgm ${img}/L2_BAND${band}-r025.png
      done
    done

  Piping the output through 
  
    gawk \
        ' BEGIN { split("",tbd); split("",tim); split("",tpx); } 
          /image [=]/ { im=$3; bd=$6; tbd[bd]=1; tim[im]=1; next;}
          /^[ ]*255[ ]/ { tpx[bd,im]=$2; next; }
          END { 
            nim = asorti(tim); nbd = asorti(tbd);
            printf "  "; 
            for (j = 1; j <= nim; j++) { printf " %8s", tim[j]; } 
            printf "\n";
            for (i = 1; i <= nbd; i++) {
              printf "%d ", tbd[i]; 
              for (j = 1; j <= nim; j++) { printf " %8d", tpx[tbd[i],tim[j]]; } 
              printf "\n";
            }
          }
        ' \
      | sort -b -k2,2n
      

             L52a     L53a     L62a     L63a
      6       200      323     2435      847
      4     18319      502    26741      405
      5     27355      221    44369     4833
      7     33428      243    50107     2238
      3     51161      218    71301     1456
      2     51715       40    70923      574
      1     89691       61   109559     1575

  Selecting "${img}/t7-th-max.pgm as the cloud mask:
  
    for img in L52a L62a L53a L63a ; do 
      band=7
      tmax="`echo ${vmax[${band}]}/255 | bc -lq`"
      echo "image = ${img} band = ${band} threshold = ${vmax[${band}]} = ${tmax}"
      convert ${img}/L2_BAND${band}-r025.png PGM:- \
        | pnmdepth 255 \
        | pamthreshold -simple -threshold="${tmax}" \
        | pamtopnm \
        | pnmdepth 255 \
        > ${img}/clouds-r100.png
      display -resize '25%' ${img}/clouds-r100.png
    done
    
COMPUTING IMAGE ALIGNMENT MATRICES

  Trying to stitch L52a and L62a. Collected lists of pairs of points,
  
    match-points-L52a-L62a.txt
    match-points-L53a-L63a.txt
    match-points-L52a-L53a.txt
    match-points-L62a-L63a.txt
  
  Used my program "image_stitch_n" to compute projective matrices
  (actually affine only):
  
    compute-stitching-matrices.sh
    
  Outputs are files "matrices/bm-matrix-{ID}-dir.txt" and 
  "matrices/bm-matrix-{ID}-inv.txt". 

    catsep matrices/bm-matrix*dir.txt \
      | sed \
          -e 's/[~][~]*//g' \
          -e 's/[\#]FILE *//g' \
          -e '/[\#] *Last/d' \
          -e '/[\#] *Created/d' \
          -e 's/^/      /g'

      matrices/bm-matrix-L52a-dir.txt
      
        1.0000000000000000e+00   6.6872239000193649e+03  -2.1337253315592534e+01
        0.0000000000000000e+00   1.0001480835738397e+00  -2.1948923557890687e-04
        0.0000000000000000e+00  -6.0687536784698182e-04   1.0012970640201990e+00
      
      matrices/bm-matrix-L53a-dir.txt
      
        1.0000000000000000e+00   5.5109616055952547e+03   5.3265051820648350e+03
        0.0000000000000000e+00   1.0000000000000000e+00   0.0000000000000000e+00
        0.0000000000000000e+00   0.0000000000000000e+00   1.0000000000000000e+00
      
      matrices/bm-matrix-L62a-dir.txt
      
        1.0000000000000000e+00   9.4557291463215370e+02  -2.9511893227843757e+01
        0.0000000000000000e+00   1.0046148412243781e+00   3.3526719684040529e-04
        0.0000000000000000e+00   4.7413532549045456e-04   1.0008857137054006e+00
      
      matrices/bm-matrix-L63a-dir.txt
      
        1.0000000000000000e+00  -2.2158489233711134e+02   5.3143815213615126e+03
        0.0000000000000000e+00   1.0046420722651916e+00   4.3084160713924291e-04
        0.0000000000000000e+00   1.2771849399536572e-03   9.9994826076389298e-01

  The root mean square and max matching errors are
  
    L52a L62a 0.7113  1.2638
    L53a L63a 1.0034  1.8136
    L52a L53a 0.9940  1.8951
    L62a L63a 0.9337  1.5817

  Checking the mismatch vectors:
  
    plot-mismatches.sh L52a L53a 300
    plot-mismatches.sh L52a L62a 300
    plot-mismatches.sh L53a L63a 300
    plot-mismatches.sh L62a L63a 300
  
  There is no systematic pattern in these vectors that would suggest a
  non-affine distortion. Therefore I will use these maps for the stitching.
  
  Mapped domain bounding boxes:
  
    img       xlo     xhi      ylo     yhi     xctr    yctr
    ----  ------- -------  ------- -------  ------- -------
    L52a   6904.7 14358.8      8.4  6670.9  10631.8  3339.6
    L53a   5770.0 13217.0   5353.5 11995.5   9493.5  8674.5
    L62a   1171.3  8628.3      5.0  6648.6   4899.8  3326.8
    L63a     45.0  7496.2   5346.0 11976.8   3770.6  8661.4

  The overall bounding box of the mapped domains is
  
    [     45.0 _  14358.8 ] × [      5.0 _  11995.5 ] center = (   7201.9   6000.3 )

  We can neatly round this up to 14400 (=6*2400) by 12000 (=5*2400) pixels.

PERFORMING IMAGE ALIGNMENT
    
  Using "pnmprojmap" to adjust the images:
  
    # bands=( 1 2 3 4 5 7 )
    bands=( 4 5 7 )
    for band in ${bands[@]} ; do 
      for img in L52a L62a L53a L63a ; do 
        cd ${tdir}/${img}
        convert L2_BAND${band}-r100.png PGM:- \
          | pnmprojmap \
              -xAxis right \
              -yAxis down \
              -iOrg 0 0 \
              -undef 0.0 \
              -matrix `cat ../matrices/bm-matrix-${img}-dir.txt | sed -e 's/[\#].*$//g'` \
              -oSize 14400 12000 \
              -oOrg 0 0 \
              -maxval 255 \
              -verbose \
              - \
          | convert PGM:- L2_BAND${band}-r100-m.png
        display -resize '5%' L2_BAND${band}-r100-m.png
      done
    done
      
CHECKING IMAGE ALIGNMENT
    
  Pairwise bounding box intersections:
    
    gawk '
      BEGIN{
        split("", Xlo); split("", Xhi); split("", Ylo); split("", Yhi); 
        Xlo["L52a"] = 6904.7; Xhi["L52a"] = 14358.8; Ylo["L52a"] =    8.4; Yhi["L52a"] =  6670.9;
        Xlo["L53a"] = 5770.0; Xhi["L53a"] = 13217.0; Ylo["L53a"] = 5353.5; Yhi["L53a"] = 11995.5;
        Xlo["L62a"] = 1171.3; Xhi["L62a"] =  8628.3; Ylo["L62a"] =    5.0; Yhi["L62a"] =  6648.6;
        Xlo["L63a"] =   45.0; Xhi["L63a"] =  7496.2; Ylo["L63a"] = 5346.0; Yhi["L63a"] = 11976.8;
      }
      //{ p1=$1; p2=$2; 
        xlo1=Xlo[p1]; xhi1=Xhi[p1]; ylo1=Ylo[p1]; yhi1=Yhi[p1];
        xlo2=Xlo[p2]; xhi2=Xhi[p2]; ylo2=Ylo[p2]; yhi2=Yhi[p2];
        xlo=int((xlo1 > xlo2 ? xlo1 : xlo2) - 1);
        xhi=int((xhi1 < xhi2 ? xhi1 : xhi2) + 1);
        ylo=int((ylo1 > ylo2 ? ylo1 : ylo2) - 1);
        yhi=int((yhi1 < yhi2 ? yhi1 : yhi2) + 1);
        nx = xhi + 1 - xlo;
        ny = yhi + 1 - ylo;
        printf "%6s%-5s %-5s  %+6d %+6d %+6d %+6d  %+6d %+6d %+6d %+6d  %05dx%05d+%05d+%05d\n", \
          "", p1, p2, xlo1, xhi1, ylo1, yhi1, xlo2, xhi2, ylo2, yhi2,  nx, ny, xlo, ylo;
      } '

      L52a  L53a    +6904 +14358     +8  +6670   +5770 +13217  +5353 +11995  06316x01320+06903+05352
      L52a  L62a    +6904 +14358     +8  +6670   +1171  +8628     +5  +6648  01727x06643+06903+00007
      L53a  L63a    +5770 +13217  +5353 +11995     +45  +7496  +5346 +11976  01729x06626+05769+05352
      L62a  L63a    +1171  +8628     +5  +6648     +45  +7496  +5346 +11976  06328x01305+01170+05345

  Let's extract those pieces and compare them:
    
    cd ${tdir}
    for band in 4 ; do 
      for pdt in \
          L52a:L53a:06316x01320+06903+05352 \
          L52a:L62a:01727x06643+06903+00007 \
          L53a:L63a:01729x06626+05769+05352 \
          L62a:L63a:06328x01305+01170+05345 \
        ; do

        tmp=( `echo ${pdt} | tr ':' ' '` )
        img1="${tmp[0]}"
        img2="${tmp[1]}"
        crop="${tmp[2]}"
        echo "comparing ${img1} and ${img2}"
        iname="L2_BAND${band}-r100-m"
        oname="L2_BAND${band}-r100-m-c${crop}"
        cname="b${band}-r100-m-c${crop}"
        for img in ${img1} ${img2} ; do 
          echo "  cropping ${img} ..."
          convert \
            ${img}/${iname}.png \
            -crop "${crop}" \
            ${img}/${oname}.pgm 
        done

        echo "  comparing ..."
        pnmxarith \
            -subtract -scale 4.0 -offset 0.5 \
            ${img1}/${oname}.pgm \
            ${img2}/${oname}.pgm \
          | pnmtopng \
          > cmp-${img1}-${img2}-${cname}.png
        display cmp-${img1}-${img2}-${cname}.png
      done
    done

  There are mismatches of up to 2 pixels at the 100% scale.
  Perhaps due to lack of perspective correction?
  Must improve image_stitch and use auto-alignment.

>>> STOPPED HERE 2012-01-04 >>>

MERGING THE IMAGES - PRELIMINARY

  Combining the two images L52a L62a:
    
    cd ${tdir}
    
    # bands=( 1 2 3 4 5 7 )
    bands=( 4 5 7 )
    for band in ${bands[@]} ; do 
      for img in L52a L62a ; do
        fname="${img}/L2_BAND${band}-r100"
        convert ${fname}-p.png ${fname}-p.pgm 
      done

      fname1="L52a/L2_BAND${band}-r100"
      fname2="L62a/L2_BAND${band}-r100"
      fname12="L52a+L62a/L2_BAND${band}-r100"
      pnmxarith \
        -maximum \
        ${fname1}-p.pgm \
        ${fname2}-p.pgm \
        | convert PGM:- ${fname12}-p.png
      rm -f ${fname1}-p.pgm
      rm -f ${fname2}-p.pgm
    done
    
    cd ${tdir}
    for band in ${bands[@]} ; do 
      fname12="L52a+L62a/L2_BAND${band}"
      convert ${fname12}-r100-p.png -resize '50%' ${fname12}-r050-p.png 
      convert ${fname12}-r100-p.png -resize '25%' ${fname12}-r025-p.png 
      convert ${fname12}-r100-p.png -resize  '5%' ${fname12}-r005-p.png 
    done

  Combining the two images L53a+L63a:
    
    cd ${tdir}
    mkdir -p L53a+L63a

    # bands=( 1 2 3 4 5 7 )
    bands=( 4 5 7 )
    for band in ${bands[@]} ; do 
      for img in L53a L63a ; do
        fname="${img}/L2_BAND${band}-r100"
        convert ${fname}-p.png ${fname}-p.pgm 
      done

      fname1="L53a/L2_BAND${band}-r100"
      fname2="L63a/L2_BAND${band}-r100"
      fname12="L53a+L63a/L2_BAND${band}-r100"
      pnmxarith \
        -maximum \
        ${fname1}-p.pgm \
        ${fname2}-p.pgm \
        | convert PGM:- ${fname12}-p.png
      rm -f ${fname1}-p.pgm
      rm -f ${fname2}-p.pgm
    done
    
    cd ${tdir}
    for band in ${bands[@]} ; do 
      fname12="L53a+L63a/L2_BAND${band}"
      convert ${fname12}-r100-p.png -resize '50%' ${fname12}-r050-p.png 
      convert ${fname12}-r100-p.png -resize '25%' ${fname12}-r025-p.png 
      convert ${fname12}-r100-p.png -resize  '5%' ${fname12}-r005-p.png 
    done
      
CHECKING THE LUT CORRECTIONS
    
    cd ${tdir}
    for band in 4 ; do 
      for uc in \
          '6703+3100' \
          '7045+3244' \
          '6862+4488' \
          '6131+5821' \
          '6621+5800' \
        ; do \
        for img in L52a L62a ; do
          fname="${img}/L2_BAND${band}-r100"
          echo "${fname}"
          convert \
              ${fname}-p.png \
              -crop "48x48+${uc}" \
              PGM:- \
            > ${fname}-p-s-"${uc}".pgm 
          cat ${fname}-p-s-"${uc}".pgm \
            | pnmnoraw \
            | sed -e '/^[#]/d' \
            | tr ' ' '\012' | sed -e '/^$/d' | tail -n +5 \
            > ${fname}-p-s-"${uc}".dat
        done
        display -filter Box -resize '400%' -contrast-stretch '1x1%' {L52a,L62a}/L2_BAND${band}-r100-p-s-"${uc}".pgm

        fname1="L52a/L2_BAND${band}-r100-p-s-"${uc}""
        fname2="L62a/L2_BAND${band}-r100-p-s-"${uc}""
        fname12="L52a+L62a/L2_BAND${band}-r100-p-s-"${uc}""
        pr --merge --omit-pagination \
          ${fname1}.dat \
          ${fname2}.dat \
          > ${fname12}.plt
        gnuplot \
          -e "set size ratio -1" \
          -e "plot \"${fname12}.plt\" using 1:2 with points pt 7 ps 0.5 lc rgb '#0044ff'" \
          -e "pause 300" &
        fit-line.gawk ${fname12}.plt > ${fname12}.parms 
        cat ${fname12}.parms
      done
    done
      

COMBINING CHANNELS

    cd ${tdir}/L52a+L62a
    
    # bands=( 1 2 3 4 5 7)
    bands=( 4 5 7)
    
    declare -a vmin
    vmin[1]=43
    vmin[2]=14
    vmin[3]=8
    vmin[4]=0
    vmin[5]=0
    vmin[6]=116
    vmin[7]=1
    
    declare -a vmax
    vmax[1]=90
    vmax[2]=50
    vmax[3]=50
    vmax[4]=118
    vmax[5]=103
    vmax[6]=144
    vmax[7]=49

    chr=5
    chg=4
    chb=7
    for band in ${chr} ${chg} ${chb} ; do
      convert L2_BAND${band}-r025-p.png PGM:- \
        | pnmdepth 255 \
        | pnmnorm \
            -bvalue=${vmin[${band}]} \
            -wvalue=${vmax[${band}]} \
        > t${band}.pgm
      display -resize '25%' t${band}.pgm
    done
    rgb3toppm t${chr}.pgm t${chg}.pgm t${chb}.pgm \
      | convert PGM:- t${chr}${chg}${chb}.png

  The 547 combination seems to give the most "natural" results, with medium-green forest and 
  light green/beige/brick farmland.