program rw2pnm C..... C This program converts "raw" binary files to pgm or ppm files. C Log compression can be used. C..... implicit real (a-h,o-z) character base*1024, filein*1024, filout*1024, & frmt*20, ccol*4, crow*4, tail(2)*4 logical input, output data tail /'.pgm', '.ppm'/ C..... C h,s,v = hue, saturation, and brightness (or "value") for pixel C ir,ig,ib = corresponding red, green, and blue values for pixel C eps = "small" number added to data to prevent underflows C xnorm = number used to normalize all input C xpnt = the current field value obtained from the input file C i,j = dummy counters C iin,iout = record pointers for the input and output files C irow,icol= number of rows and columns of data (integers) C crow,ccol= number of rows and columns of data (characters) C nframe = current frame number C ndigit = number of decimal digits in the current frame number C ndec = number of decades of compression C input,output = true if input/output file exists (logical) C base = base input and output file names C filein,filout = complete input and output file names C frmt = internal file used to construct file names C iside = indicates two-sided or one-sided data C icolor = indicates color or grayscale C tail = contains possible extensions for output file names C xtmp = temporary variable for holding result C..... write(6,*)'Enter the base name of the raw files.' read(5,'(a)')base C.....Determine the location of last character in the base name. ilen = index(base,' ') - 1 10 write(6,*)'Enter the normalization.' read(5,*)xnorm if (xnorm .le. 0.) goto 10 20 write(6,*)'Enter the number of decades (0 = linear scale).' read(5,*)ndec if (ndec .lt. 0) goto 20 30 write(6,*)'Do you want color or grayscale output? (1 = ', & 'grayscale; 2 = color)' read(5,*)icolor if (icolor .eq. 2) then 40 write(6,*)'Do you want to obtain "one-sided" output by', & ' taking the absolute value' write(6,*)'of the input data? (1 = one-sided; 2 = two-sided)' read(5,*)iside if (iside .ne. 1 .and. iside .ne. 2) goto 40 elseif (icolor .ne. 1) then goto 30 endif nframe = 0 1000 nframe = nframe + 1 C........Calculate the number of digits in the current frame number. ndigit = log10(real(nframe)+0.01) + 1 C........Construct the proper format field for the input file name. write(frmt,'(a2,i2,a5,i9,a1)')'(a',ilen,',a1,i',ndigit,')' C........Put this together to obtain the input file name. write(filein,frmt)base,'.',nframe C........Construct the output file name. filout = filein filout(ilen+ndigit+2:ilen+ndigit+6) = tail(icolor) C........Check if input file is available. Terminate if not. inquire(file=filein,exist=input) if (.not. input) then write(6,*)'Input not found: ',filein(1:index(filein,' ')) stop endif C........If output file already exists, skip to next file. inquire(file=filout,exist=output) if (output) then write(6,*)'Skipping:',filout(1:index(filout,' ')) goto 1000 endif write(6,*)'Processing: ',filein(1:index(filein,' ')) C........Open the input file and read the size of the image. open(55,file=filein,form='unformatted',status='old', & access='direct',recl=4) read(55,rec=1)xtmp icol = int(xtmp)+10 read(55,rec=2)xtmp irow = int(xtmp) write(ccol,'(i4)')icol write(crow,'(i4)')irow C........Open the output file and write the header information. open(66,file=filout,form='unformatted',status='new', & access='direct',recl=1) write(66,rec=1)'P' write(66,rec=2)char(52+icolor) write(66,rec=3)char(10) write(66,rec=4)ccol(1:1) write(66,rec=5)ccol(2:2) write(66,rec=6)ccol(3:3) write(66,rec=7)ccol(4:4) write(66,rec=8)' ' write(66,rec=9)crow(1:1) write(66,rec=10)crow(2:2) write(66,rec=11)crow(3:3) write(66,rec=12)crow(4:4) write(66,rec=13)char(10) write(66,rec=14)'2' write(66,rec=15)'5' write(66,rec=16)'5' write(66,rec=17)char(10) C........Initialize record pointers iout=18 iin =3 C........Calculate "epsilon" that is added to values if compressing. if (ndec.ne.0) eps = 10.0**(-ndec)*xnorm/(1.0-10.0**(-ndec)) C........Loop through all pixels in the image. do 200 j=1,irow do 100 i=1,icol if (i .le. 10) then C..............Draw the color map if within the first 10 columns. if (icolor .eq. 1 .or. iside .eq. 1) then xpnt = real(irow-j)/real(irow-1)*xnorm else xpnt = (-2./real(irow-1)*real(j) & + real(irow+1)/real(irow-1))*xnorm endif else C..............If not in first 10 columns, get data from file. read(55,rec=iin,err=2000)xpnt iin = iin+1 endif C..............If the point is out of range, whack it down to size. if (abs(xpnt) .gt. xnorm) xpnt = xpnt/abs(xpnt)*xnorm C..............Use log scaling if requested. if (ndec .ne. 0) then xtmp = log10((abs(xpnt)+eps)/(xnorm+eps))/ & real(ndec)+1. else xtmp = abs(xpnt/xnorm) endif if (icolor .eq. 1) then C.................Grayscale output. write(66,rec=iout)char(int(255*xtmp)) iout = iout+1 else if (iside .eq. 1) then C.................Use color map #1 for one-sided data. call cmap1(h,s,v,xtmp) else C.................Reconstruct sign & use color map #2 for two-sided. if (xpnt .lt. 0.) xtmp = -xtmp call cmap2(h,s,v,xtmp) endif C.................Convert (h,s,v) to (r,g,b) and write output. call hsvrgb(h,s,v,ir,ig,ib) write(66,rec=iout)char(ir) write(66,rec=iout+1)char(ig) write(66,rec=iout+2)char(ib) iout = iout+3 endif 100 continue 200 continue 2000 close(55) close(66) goto 1000 end C********************************************************************** C.....Color mapping appropriate for one-sided data. subroutine cmap1(h,s,v,xpnt) real h, s, v, xpnt, pi parameter (pi=3.141593) if (xpnt .gt. 0.8) then h = 330. + 30*(1.-xpnt)/0.2 s = (1.-xpnt)/0.2 v = 1.0 elseif (xpnt .gt. 0.2) then t = 2.0*pi*((xpnt-0.2)/.6) h = 240.0*(0.8-xpnt)/0.6 + 5.0*sin(t) s = 1.0 v = 1.0 else h = 270.0*(.2-xpnt)+240.0 s = 1.0 v = xpnt/0.2 endif return end C********************************************************************** C.....Color mapping appropriate for two-sided data. Zero = black. subroutine cmap2(h,s,v,xpnt) real h, s, v, xpnt if (xpnt .gt. 0.8) then h = 30. - 30.*(xpnt-0.8)*(xpnt-0.8)/0.04 s = (1.0-xpnt)/0.2 v = 1.0 else if (xpnt .gt. .3) then h = 135.0 - 105.*(xpnt-.3)/.5 s = 1.0 v = 1.0 else if (xpnt .gt. 0.) then h = 180.0 - 45.*xpnt/.3 s = 1.0 v = xpnt/.3 else if (xpnt .gt. -.3) then h = (180. - 45.*xpnt/.3) s = 1.0 v = abs(xpnt)/.3 else h = 225. - 135.*(xpnt+.3)/.7 s = 1.0 v = 1.0 endif return end C********************************************************************** C.....Subroutine to convert (h,s,v) values to (r,g,b) values. subroutine hsvrgb(h,s,v,ir,ig,ib) implicit real (a-h,o-z) real pqvt(4) integer irpnt(6), igpnt(6), ibpnt(6) data irpnt/3,2,1,1,4,3/,igpnt/4,3,3,2,1,1/,ibpnt/1,1,4,3,3,2/ C..... C input: 0 <= h <= 360 (hue) C 0 <= v <= 1.0 (value or brightness) C 0 <= s <= 1.0 (saturation) C output: 0 <= ir,ig,ib <= maxval C..... if (s .eq. 0.0) then ir = 255.0*v ig = ir ib = ir return endif if (h .eq. 360.0) h=0.0 h = h / 60.0 i = int(h) f = h-i pqvt(3) = 255.0*v pqvt(1) = pqvt(3)*(1.0-s) pqvt(2) = pqvt(3)*(1.0-s*f) pqvt(4) = pqvt(3)*(1.0-s*(1.0-f)) i = i + 1 ir = pqvt(irpnt(i)) ig = pqvt(igpnt(i)) ib = pqvt(ibpnt(i)) return end C**********************************************************************