! mozi-koodo: EUC
! 05/11/17版
! 二値化またはグレスケール化した画像の中の黒い部分の面積比、
! 散らばり具合を求める。
! スキャナーで二値化(線画とか?)画像として保存したセロテープ試験
! の画像(jpgとかbmp)をgThumbとかで開いて規定の大きさに切り取り、
! ひとまずpng形式で保存。それをGIMPで開き直してアスキー形式の
! ppmで保存する(pgmやpbmの方がいいだろうが、GIMPで簡単に作れる
! という意味でppmを用いておく)。そして、
!   senyuu < hoge.ppm
! のように入力する。
! ----------------------------------------
! 但し、プログラムのデータ入力部分は、
! ppmファイルの各行のデータの個数がどういうふうになっているかに
! 応じて、書き換えなければならない。
! 以下は、GIMPでアスキー形式のppmフォーマットで保存したファイル
! (RGBデータが1行に1個ずつ書き込まれている)からの読み込みに対応させた例
!
        program tirabari
       implicit real*8 (a-h, o-z)
       dimension rgb3(2880000),r(2880000),g(2880000),b(2880000),
     & zy(2880000),zi(2880000),zq(2880000)
       character p3*2,mozi*44
          read(*,100) p3
100      format(a2)
!           write(*,*) ' マジックナンバー：',p3
          read(*,*) mozi
          read(*,*) m,n
!          write(*,*) '     幅：',m
!          write(*,*) '   高さ：',n
           mn=m*n
          read(*,*) max
!          write(*,*) ' 最大値：',max
          do i=1,mn
           read(*,*) rgb3((i-1)*3+1)
           read(*,*) rgb3((i-1)*3+2)
           read(*,*) rgb3((i-1)*3+3)
          end do
!
        do i=1,mn
          r(i)=rgb3((i-1)*3+1)
          g(i)=rgb3((i-1)*3+2)
          b(i)=rgb3((i-1)*3+3)
       end do
           kuro=0
!
             nr1=r(1)
          do i=1,mn
             nr=r(i)
! 黒い部分の個数を数える(グレースケールだからR値だけ数える)
! スキャナーとかで2値化しても、完全にRGB=(0,0,0)にはならないので、
! 閾値を10以下とかに設定しておく。
             nsikii=10
           if(nr.le.nsikii) kuro=kuro+1
          end do
!
!          print*,'   閾値： ',nsikii
!
             ntigau=0
          do j=1,n
             mj1=m*(j-1)+1
             nr1=r(mj1)
          do i=mj1,mj1+m-1
             nr=r(i)
             if(nr.ne.nr1) ntigau=ntigau+1
             nr1=nr
          end do
          end do
!
!
!           write(*,*) '         全ピクセル数：',mn
!           write(*,*) '         黒ピクセル数：',kuro
!           write(*,*) ' 左右色違いピクセル数：',ntigau
            nhasi=moku+nhoka
!       write(*,200) '     黒/全体(サビ面積比)：',real(kuro)/real(mn)
!       write(*,200) ' 色違い/全体(散らばり度)：',real(ntigau)/real(mn)
             hi=real(kuro)/real(mn)
             ti=real(ntigau)/real(mn)
200        format(a, f5.3)
!
!
! RGB の平均を YIQ に変換
        do i=1,mn
         zy(i)=0.31d0 * r(i) + 0.59d0 * g(i) + 0.11d0 * b(i)
!        zi(i)=0.60d0 * r(i) - 0.28d0 * g(i) - 0.32d0 * b(i)
!        zq(i)=0.21d0 * r(i) - 0.52d0 * g(i) + 0.31d0 * b(i)
        end do
          call heikin(mn,zy,ay,sy)
!         call heikin(mn,zi,ai,si)
!         call heikin(mn,zq,aq,sq)
!        write(*,*) 'Y値の平均=',ay,'（標準偏差：',sy,'）'
!       write(*,*) 'I値の平均=',ai,'（標準偏差：',si,'）'
!       write(*,*) 'Q値の平均=',aq,'（標準偏差：',sq,'）'
!        write(*,*) '------------------------------------------'
        x=real(max)
!        write(*,*) '最大値1に正規化'
!        write(*,*) 'Y値の平均=',ay/x,'（標準偏差：',sy/x,'）'
        yh=ay/x
        hy=sy/x
!       write(*,*) 'I値の平均=',ai/x,'（標準偏差：',si/x,'）'
!       write(*,*) 'Q値の平均=',aq/x,'（標準偏差：',sq/x,'）'
!
        write(*,300) hi,ti,yh,hy
300     format(f6.3,',',f6.3,',',f6.3,',',f6.3)
          end
!
!
          subroutine heikin(mn,x,a,s)
       implicit real*8 (a-h, o-z)
       dimension x(mn)
         a=0.d0
         s=0.d0
        do i=1,mn
         a=a+x(i)
        end do
         a=a/float(mn)
        do i=1,mn
         s=s+ ( x(i)-a )**2
        end do
         s=sqrt(s/mn)
        return
        end
!
