• View
• Edit
• History
• Print

Annex

Code source de l’algorithme pour obtenir les 50 pixels les plus sombres des images obtenues du radiomètre

c programme convertit un fichier fits en pgm c c c Copyright (C) 2011 Martin Aube c c This program is free software: you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation, either version 3 of the License, or c (at your option) any later version. c c This program is distributed in the hope that it will be useful, c but WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the c GNU General Public License for more details. c c You should have received a copy of the GNU General Public License c along with this program. If not, see <http://www.gnu.org/licenses/>. c c Contact: martin.aube@cegepsherbrooke.qc.ca c

        real lambda,out(400000),vect(400000),quad,moy
        real x(400000),XSMOM4,y(15000)
        integer image(765,510),imgout(765,510),moya,moyo
        integer histo(65535),h,s,l,m,fenetre,nbon,n,Im,nbh
        integer i,j,pos,longueur,nbdata,minimum,k,IWRITE
        do i=1,765
          do j=1,510
              image(i,j)=0
              imgout(i,j)=0
          enddo
        enddo
        n=0
        open(unit=1,file='imagetxt.tmp',status='unknown')
           longueur=764+766*510
           read(1,*) (vect(i),i=1,longueur)
        close(unit=1)
        minimum=65536
        moya=0

        do h=1,65535
            histo(h)=0
        enddo
        print*,'remplissage de la matrice et de l histogramme'
        k=0
        do i=1,765
           do j=1,510
              k=k+1
              pos=i+766*j

c out(n)=out(n)+vect(pos)

              image(i,j)=vect(pos)
              x(k)=vect(pos)
              if (k.le.15000) y(k)=vect(pos)
              if ((image(i,j).lt.minimum).and.(image(i,j).ne.0)) then
                minimum=image(i,j)
              endif
              moya=moya+image(i,j)
              histo(image(i,j))=histo(image(i,j))+1
           enddo
        enddo

c calcul du kurtosis ( coefficient d applatissement c IWRITE=1 c call STMOM4(x,k,IWRITE,XSMOM4) c print*,'kurtosis=',XSMOM4 c calcul du spectre spatial (fourier) c call FOURIE(y,15000) c ecriture du premier histogramme

       open(unit=4,file='histo.initial',status='unknown')
          do h=1,65535
            write(4,*) h,histo(h)
          enddo
       close(unit=4)

c nbon est le nombre de donnees a moyenner

      fenetre=10
      do h=1,65535
         if (histo(h).gt.0) then
            s=0
            do l=h,h+fenetre
                s=s+histo(l)
            enddo
            if (s.eq.histo(h)) then
              do m=h-100,h
                histo(m)=0
              enddo
            endif
         endif
      enddo
       do nbon=50,50
       n=0

c nbon=500

       h=0
       nbh=0

c nbh est l addition des valeures poderees par leur frequence

       do while (n.le.nbon) 
          h=h+1    
          if (histo(h).ne.0) then
                  n=n+histo(h)
                  nbh=nbh+histo(h)*h
                do i=1,765
                  do j=1,510
                    if (image(i,j).eq.h) then
                        imgout(i,j)=image(i,j)

c print*,h,i,j

                    endif
                  enddo
                enddo
           endif
        enddo
        moy=real(nbh)/real(n)
        print*,moy, nbon



          enddo

c ecriture du nouveau histogramme

       open(unit=3,file='histo.final',status='unknown')
          do h=1,65535
            write(3,*) h,histo(h)
          enddo
       close(unit=3)
        nbdata=765*510











        print*,'ecriture des pgm'
	open(unit=2,file='image.pgm',status='unknown')
          write(2,100)

c write(2,200) int(1.5*moy)

          write(2,200) 



          write(2,*) ((image(i,j),i=1,765),j=1,510)
         close(unit=2)
	open(unit=1,file='image1.pgm',status='unknown')
          write(1,100)

c write(1,200) int(1.*moy)

          write(1,200) 


          write(1,*) ((imgout(i,j),i=1,765),j=1,510)
         close(unit=1)

 100     format('P2')

c 200 format('765 510 ',I4 )

 200     format('765 510 65535' )
	stop
	end