!------------------------------------------------------------------------ ! PRODUCES A BMP IMAGE FROM A RADMC-3D image.out FILE ! ! The pixout() subroutine and all dependent subroutines were written by ! Keiji Hayashi, and this software can be found at: ! ! http://sun.stanford.edu/~keiji/pixelfrt.html ! ! Thanks, Keiji, for making this available!! !------------------------------------------------------------------------ module image_to_bmp double precision, allocatable :: image(:,:) ! RADMC-3D image data array character*1, allocatable :: bitmap(:,:,:) ! RGB data array integer :: nx,ny integer :: bnx,bny integer :: factor contains !------------------------------------------------------------------------ ! THE MAIN ROUTINE FOR THIS PROGRAM !------------------------------------------------------------------------ subroutine convert_image_to_bmp() implicit none double precision :: mn,mx integer :: ilg,ict integer :: ct_r(0:255),ct_g(0:255),ct_b(0:255) logical :: lg,ct ! ! Ask some things ! write(*,*) 'Minimum intensity?' read(*,*) mn write(*,*) 'Maximum intensity?' read(*,*) mx write(*,*) 'Log? (1=yes, 0=no)' read(*,*) ilg if(ilg.eq.0) then lg=.false. else lg=.true. endif write(*,*) 'Size magnification factor?' read(*,*) factor write(*,*) 'Use color table? (1=yes: read from ct.inp, 0=no: use greyscale)' read(*,*) ict if(ict.eq.0) then ct=.false. else ct=.true. endif ! ! If required, read the color table ! if(ct) then open(unit=1,file='ct.inp') read(1,*) ct_r read(1,*) ct_g read(1,*) ct_b close(1) endif ! ! Read the image ! call read_radmc3d_image() ! ! Allocate the bitmap ! bnx = nx*factor bny = ny*factor allocate(bitmap(3,bnx,bny)) ! ! Create the bitmap ! if(ct) then call convert_monoimage_to_bitmap(nx,ny,bnx,bny,image,mn,mx,lg,bitmap, & ct_r,ct_g,ct_b) else call convert_monoimage_to_bitmap(nx,ny,bnx,bny,image,mn,mx,lg,bitmap) endif ! ! Write the bitmap to BMP ! call pixout(bnx,bny,bitmap) ! ! Cleanup ! if(allocated(image)) deallocate(image) if(allocated(bitmap)) deallocate(bitmap) end subroutine convert_image_to_bmp !------------------------------------------------------------------------ ! READ RADMC-3D IMAGE !------------------------------------------------------------------------ subroutine read_radmc3d_image() implicit none integer :: iformat,i,j double precision :: dum(3) open(unit=1,file='image.out') read(1,*) iformat if(iformat.ne.1) stop "Unknown RADMC-3D image format" read(1,*) nx,ny read(1,*) iformat read(1,*) dum(1:3) allocate(image(nx,ny)) read(1,*) ((image(i,j),i=1,nx),j=1,ny) close(1) end subroutine read_radmc3d_image !------------------------------------------------------------------------ ! CREATE BITMAP FROM IMAGE !------------------------------------------------------------------------ subroutine convert_monoimage_to_bitmap(nx,ny,bnx,bny,image,mn,mx,lg,rgb, & ct_r,ct_g,ct_b) implicit none integer :: nx,ny,bnx,bny integer, optional :: ct_r(0:255),ct_g(0:255),ct_b(0:255) double precision :: image(nx,ny) double precision :: mn,mx logical :: lg character*1 rgb(3,bnx,bny) ! RGB data array double precision :: scaledint integer :: i,j,bi,bj,ired,igre,iblu,icol ! do bj = 1, bny do bi = 1, bnx ! ! Find the i,j ! i = bi / factor j = bj / factor ! ! Create the scaled intensity ! if(lg) then scaledint = ( log(image(i,j)+1d-90) - log(mn) ) / & ( log(mx) - log(mn) ) else scaledint = ( image(i,j) - mn ) / & ( mx - mn ) endif if(scaledint.gt.1.d0) scaledint = 1.d0 if(scaledint.lt.0.d0) scaledint = 0.d0 ! ! Create the three colors ! if(present(ct_r)) then ! ! Color table ! icol = int(scaledint * 255.0D+00) if (icol .GT. 255) icol = 255 if (icol .LT. 0) icol = 0 ired = ct_r(icol) igre = ct_g(icol) iblu = ct_b(icol) else ! ! Greyscale ! ired = int(scaledint * 255.0D+00) igre = int(scaledint * 255.0D+00) iblu = int(scaledint * 255.0D+00) endif ! ! Check range ! if (ired .GT. 255) ired = 255 if (igre .GT. 255) igre = 255 if (iblu .GT. 255) iblu = 255 if (ired .LT. 0) ired = 0 if (igre .LT. 0) igre = 0 if (iblu .LT. 0) iblu = 0 ! ! Put into bitmap ! rgb(1,bi,bj) = char(ired) rgb(2,bi,bj) = char(igre) rgb(3,bi,bj) = char(iblu) enddo enddo end subroutine convert_monoimage_to_bitmap !======================================================================== ! THE SUBROUTINES BELOW WERE WRITTEN BY KEIJI HAYASHI ! ! http://sun.stanford.edu/~keiji/pixelfrt.html ! ! Converted into f90 by C. Dullemond 04.06.2012 ! !======================================================================== !------------------------------------------------------------------------ ! ! Notes ! o With a parameter ipixout set at 1, 2 or others, ! this subroutine will generate PPM-P6(binary), PPM-P3(text), ! or BMP(24bit depth without color table). ! ! o Some parts follow DEC-FORTRAN convention that had been defacto-standard long ago. ! Some compilers today may not accept if "ipixout" is set other than 2. ! ! o g77 (ver. 3.3.3) works for all three choices. ! o Intel compiler (ver. 9 or later) works for all three choices. ! ! -------------------------------------------- subroutine pixout(ihpixf,jvpixf,rgb) implicit none ! interface arg. integer ihpixf, jvpixf !parameter(ihpixf = 128, jvpixf = 128) ! pixel size, eacg must be multiple of 4, if BMP is chosen as output format. character*1 rgb(3,ihpixf,jvpixf) ! RGB data array ! local character*12 fnameout integer i, j, k integer itmp, icnt character*14 frmtstr character*54 headmsw character*4 byt4 character*2 byt2 ! choices integer ipixout parameter(ipixout = 3) ! 1 / 2 / other= PPM6, PPM3, BMP(24bit) if (ipixout .EQ. 1) then ! ! PPM P6 ! fnameout = 'image.ppm' open(unit=2,file=fnameout,status='unknown') write(*,*) 'Now writing PPM (P6) file : ', fnameout ! ! header ! write(2,'(''P6'', 2(1x,i4),'' 255 '',$)') & ! some compiler may not accept this line. ihpixf, jvpixf ! ! image data ! itmp = ihpixf * jvpixf * 3 write(frmtstr,'(''('',i8.8,''A,$)'')') itmp ! make output "format" write(2,fmt=frmtstr) & ! some compiler may not accept this line. (((rgb(k,i,j),k=1,3),i=1,ihpixf),j=jvpixf,1,-1) ! here, j (vertical address) runs from top to bottom. close(2) ! else if (ipixout .EQ. 2) then ! ! PPM P3 ! rather "safer" choice for many Fortran compiler(s). ! fnameout = 'image.ppm' open(unit=2,file=fnameout,status='unknown') write(*,*) 'Now writing PPM (P3) file : ', fnameout ! ! header ! write(2,'(A)') 'P3' write(2,'(2(1x,i4),'' 255 '')') ihpixf, jvpixf icnt = 0 ! ! image data ! do j = jvpixf, 1, -1 ! here, j (vertical address) runs from top to bottom. do i = 1, ihpixf, 1 do k = 1, 3 itmp = ichar(rgb(k,i,j)) icnt = icnt + 4 if (icnt .LT. 60) then write(2,fmt='(1x,i3,$)') itmp ! "$" is not standard. else write(2,fmt='(1x,i3)') itmp icnt = 0 endif enddo enddo enddo write(2,'(A)') ' ' close(2) ! else ! ! BMP (24bit depth)... this part works only when width is multiple of 4. ! itmp = mod(ihpixf, 4) if (itmp .NE. 0) then write(*,*) 'width must be multiple of 4' stop endif ! fnameout = 'image.bmp' open(unit=2,file=fnameout,status='unknown') write(*,*) 'Now writing BMP(24bit) file : ', fnameout ! ! header 1 (file header ; 1--14 byte) ! headmsw( 1: 2) = 'BM' ! declaring this is BMP file itmp = 54 + ihpixf * jvpixf * 3 ! total file size = header + data call num2bit4(itmp,byt4) headmsw( 3: 6) = byt4(1:4) itmp = 0 ! may be 0 call num2bit2(itmp,byt2) headmsw( 7: 8) = byt2(1:2) itmp = 0 ! may be 0 call num2bit2(itmp,byt2) headmsw( 9:10) = byt2(1:2) itmp = 54 ! must be 54 : total length of header call num2bit4(itmp,byt4) headmsw(11:14) = byt4(1:4) ! ! header 2 (bit-map header ; 13--54 byte) ! itmp = 40 ! must be 40 : length of bit-map header call num2bit4(itmp,byt4) headmsw(15:18) = byt4(1:4) itmp = ihpixf ! width call num2bit4(itmp,byt4) headmsw(19:22) = byt4(1:4) itmp = jvpixf ! height call num2bit4(itmp,byt4) headmsw(23:26) = byt4(1:4) itmp = 1 ! must be 1 call num2bit2(itmp,byt2) headmsw(27:28) = byt2(1:2) itmp = 24 ! must be 24 : color depth in bit. call num2bit2(itmp,byt2) headmsw(29:30) = byt2(1:2) itmp = 0 ! may be 0 : compression method index call num2bit4(itmp,byt4) headmsw(31:34) = byt4(1:4) itmp = 0 ! may be 0 : file size if compressed call num2bit4(itmp,byt4) headmsw(35:38) = byt4(1:4) itmp = 0 ! arbit. : pixel per meter, horizontal call num2bit4(itmp,byt4) headmsw(39:42) = byt4(1:4) itmp = 0 ! arbit. : pixel per meter, vertical call num2bit4(itmp,byt4) headmsw(43:46) = byt4(1:4) itmp = 0 ! may be 0 here : num. of color used call num2bit4(itmp,byt4) headmsw(47:50) = byt4(1:4) itmp = 0 ! may be 0 here : num. of important color call num2bit4(itmp,byt4) headmsw(51:54) = byt4(1:4) ! ! writing header part ! write(2,'(a54,$)') headmsw(1:54) ! ! image data ! itmp = ihpixf * jvpixf * 3 write(frmtstr,'(''('',i8.8,''A,$)'')') itmp write(2,fmt=frmtstr) & (((rgb(k,i,j),k=3,1,-1),i=1,ihpixf),j=1,jvpixf) ! writing in BGR order, not RGB. close(2) endif return end subroutine pixout !-------------------------------------- ! convert integer values to 4 8-bit characters !-------------------------------------- subroutine num2bit4(inum,byt4) implicit none integer inum character*4 byt4 integer itmp1, itmp2 itmp1 = inum itmp2 = itmp1 / 256**3 byt4(4:4) = char(itmp2) itmp1 =-itmp2 * 256**3 +itmp1 itmp2 = itmp1 / 256**2 byt4(3:3) = char(itmp2) itmp1 =-itmp2 * 256**2 +itmp1 itmp2 = itmp1 / 256 byt4(2:2) = char(itmp2) itmp1 =-itmp2 * 256 +itmp1 byt4(1:1) = char(itmp1) return end subroutine num2bit4 ! -------------------------------------- ! convert integer values to 2 8-bit characters ! -------------------------------------- subroutine num2bit2(inum,byt2) implicit none integer inum character*2 byt2 integer itmp1, itmp2 itmp1 = inum itmp2 = itmp1 / 256 byt2(2:2) = char(itmp2) itmp1 =-itmp2 * 256 + itmp1 byt2(1:1) = char(itmp1) return end subroutine num2bit2 end module image_to_bmp !======================================================================== ! THE MAIN PROGRAM !======================================================================== program main use image_to_bmp call convert_image_to_bmp() end program main