c 03/12/14版
c アスキーフォーマットのppm画像ファイルの R値、G値、B値それぞれの
c 平均と標準偏差を求める。また、
c ppm画像データの各行の輝度信号(Y=0.31R+0.59G+0.11B)の横方向の平均を
c 1列に並べたデータ（yokosima.d）と、
c ppm画像データの各列の輝度信号(Y=0.31R+0.59G+0.11B)の縦方向の平均を
c 1列に並べたデータ（tatesima.d）を出力する。ppm.fがやるのはここまで。
c ------------------意図----------------
c   貴家仁志『よくわかるディジタル画像処理』(CQ出版社)
c   付属の1次元FFTプログラム（fftvec）に、
c yokosima.d とtatesima.d を
c 入力して、画像ファイルの横縞成分の振幅スペクトルyokospe.dと
c 縦縞成分の振幅スペクトルtatespe.dを求める。
c それらをspe2dg.fに入力すると、パワーと周波数で両対数プロット
c した場合の横縞パワーと縦縞パワーの回帰直線の
c 傾きの平均、標準偏差の平均、傾きの差などを
c 定量化のための数値指標として出力する。
c 例えば、gazou01.ppm という画像データに対しては、
c ppm < gazou01.ppm > rgb.01
c fftvec yokosima.d > yokospe.d
c fftvec tatesima.d > tatespe.d
c spe2dg > spe.01
c みたいな手順で実行する。
c ----------------------------------------
c だけど、データ入力は、
c ppmファイルの各行のデータの個数がどういうふうになっているかに
c 応じて、書き換えなければならない。
c 以下は、
c http://www.vector.co.jp/soft/win95/art/se165334.html
c ここのフリーソフト「疾風 -Tokikaze-」（作者：ＰＡＰＹさん）
c で、256×256ピクセルで出力された ppm ファイルに対応した例。
c
        program ppm
       implicit real*8 (a-h, o-z)
       dimension rgb3(786432),r(262144),g(262144),b(262144),
     & zy(262144),zi(262144),zq(262144)
       character p3*2,mozi*35
          read(*,100) p3
          read(*,110) mozi
100      format(a2)
110      format(a35)
           write(*,*) 'マジックナンバー：',p3
           write(*,*) 'コメント行：',mozi
          read(*,*) m,n
          write(*,*) '幅：',m
          write(*,*) '高さ：',n
           mn=m*n
          read(*,*) max
          write(*,*) '最大値：',max
          read(*,*) (rgb3(i),i=1,18)
          n1=18
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,15)
          n1=n1+15
c          write(*,*) 'n1=',n1
          mn3=mn*3
          n2=mn3-n1-9
          n2=n2/69
c          write(*,*) 'n2=',n2
        do j=1,n2
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,18)
          n1=n1+18
          read(*,*) (rgb3(n1+i),i=1,15)
          n1=n1+15
        end do
          read(*,*) (rgb3(n1+i),i=1,9)
c          write(*,*) 'nsaigo=',n1+9
c  以上がデータ入力
c
        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
          call heikin(mn,r,ar,sr)
          call heikin(mn,g,ag,sg)
          call heikin(mn,b,ab,sb)
        write(*,*) 'R値の平均=',ar,'（標準偏差：',sr,'）'
        write(*,*) 'G値の平均=',ag,'（標準偏差：',sg,'）'
        write(*,*) 'B値の平均=',ab,'（標準偏差：',sb,'）'
           x=float(max)
        write(*,*) '------------------------------------------'
        write(*,*) '最大値1に正規化'
        write(*,*) 'R値の平均=',ar/x,'（標準偏差：',sr/x,'）'
        write(*,*) 'G値の平均=',ag/x,'（標準偏差：',sg/x,'）'
        write(*,*) 'B値の平均=',ab/x,'（標準偏差：',sb/x,'）'
        write(*,*) '------------------------------------------'
c
c 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(*,*) '------------------------------------------'
        write(*,*) '最大値1に正規化'
        write(*,*) 'Y値の平均=',ay/x,'（標準偏差：',sy/x,'）'
        write(*,*) 'I値の平均=',ai/x,'（標準偏差：',si/x,'）'
        write(*,*) 'Q値の平均=',aq/x,'（標準偏差：',sq/x,'）'
c
       open(7,file='yokosima.d')
          write(7,*) m
        do i=1,n
          im=(i-1)*m
          yoko=0.d0
         do j=im+1,im+m
          yoko=yoko+0.31d0*r(j) +0.59d0*g(j) +0.11d0*b(j)
         end do
          yoko=yoko/float(max)
          write(7,200) yoko
c          write(7,210) i,yoko
        end do
       close(7)
c
       open(7,file='tatesima.d')
          write(7,*) m
        do j=1,m
          mnj=m*(n-1)+j
          tate=0.d0
         do i=j,mnj,m
          tate=tate+0.31d0*r(i) +0.59d0*g(i) +0.11d0*b(i)
         end do
          tate=tate/float(max)
          write(7,200) tate
c          write(7,210) j,tate
        end do
       close(7)
c
200     format(1pd13.5)
210     format(i4,1pd13.5)
c
          end
c
          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
c
