The main routine used to convert the raw files to PGM or PPM files is rw2pnm. (As mention previously, it is important to note that the C version of this code, which is available from the Web site, is much simpler and significantly faster than this Fortran code.) This program relies upon the subroutines cmap1 and cmap2 to generate the one-sided and two-sided color maps, respectively. The subroutine hsvrgb is used to convert from values in HSV space to RGB values. The complete program listing for rw2pnm and its subroutines follows:
program rw2pnm 1
C..... 2
C This program converts "raw" binary files to pgm or ppm files. 3
C Log compression can be used. 4
C..... 5
implicit real (a-h,o-z) 6
character base*1024, filein*1024, filout*1024, 7
& frmt*20, ccol*4, crow*4, tail(2)*4 8
logical input, output 9
data tail /'.pgm', '.ppm'/ 10
C..... 11
C h,s,v = hue, saturation, and brightness (or "value") for pixel 12
C ir,ig,ib = corresponding red, green, and blue values for pixel 13
C eps = "small" number added to data to prevent underflows 14
C xnorm = number used to normalize all input 15
C xpnt = the current field value obtained from the input file 16
C i,j = dummy counters 17
C iin,iout = record pointers for the input and output files 18
C irow,icol= number of rows and columns of data (integers) 19
C crow,ccol= number of rows and columns of data (characters) 20
C nframe = current frame number 21
C ndigit = number of decimal digits in the current frame number 22
C ndec = number of decades of compression 23
C input,output = true if input/output file exists (logical) 24
C base = base input and output file names 25
C filein,filout = complete input and output file names 26
C frmt = internal file used to construct file names 27
C iside = indicates two-sided or one-sided data 28
C icolor = indicates color or grayscale 29
C tail = contains possible extensions for output file names 30
C xtmp = temporary variable for holding result 31
C..... 32
write(6,*)'Enter the base name of the raw files.' 33
read(5,'(a)')base 34
C.....Determine the location of last character in the base name. 35
ilen = index(base,' ') - 1 36
10 write(6,*)'Enter the normalization.' 37
read(5,*)xnorm 38
if (xnorm .le. 0.) goto 10 39
20 write(6,*)'Enter the number of decades (0 = linear scale).' 40
read(5,*)ndec 41
if (ndec .lt. 0) goto 20 42
30 write(6,*)'Do you want color or grayscale output? (1 = ', 43
& 'grayscale; 2 = color)' 44
read(5,*)icolor 45
if (icolor .eq. 2) then 46
40 write(6,*)'Do you want to obtain "one-sided" output by', 47
& ' taking the absolute value' 48
write(6,*)'of the input data? (1 = one-sided; 2 = two-sided)' 49
read(5,*)iside 50
if (iside .ne. 1 .and. iside .ne. 2) goto 40 51
elseif (icolor .ne. 1) then 52
goto 30 53
endif 54
nframe = 0 55
1000 nframe = nframe + 1 56
C........Calculate the number of digits in the current frame number. 57
ndigit = log10(real(nframe)+0.01) + 1 58
C........Construct the proper format field for the input file name. 59
write(frmt,'(a2,i2,a5,i9,a1)')'(a',ilen,',a1,i',ndigit,')' 60
C........Put this together to obtain the input file name. 61
write(filein,frmt)base,'.',nframe 62
C........Construct the output file name. 63
filout = filein 64
filout(ilen+ndigit+2:ilen+ndigit+6) = tail(icolor) 65
C........Check if input file is available. Terminate if not. 66
inquire(file=filein,exist=input) 67
if (.not. input) then 68
write(6,*)'Input not found: ',filein(1:index(filein,' ')) 69
stop 70
endif 71
C........If output file already exists, skip to next file. 72
inquire(file=filout,exist=output) 73
if (output) then 74
write(6,*)'Skipping:',filout(1:index(filout,' ')) 75
goto 1000 76
endif 77
write(6,*)'Processing: ',filein(1:index(filein,' ')) 78
C........Open the input file and read the size of the image. 79
open(55,file=filein,form='unformatted',status='old', 80
& access='direct',recl=4) 81
read(55,rec=1)xtmp 82
icol = int(xtmp)+10 83
read(55,rec=2)xtmp 84
irow = int(xtmp) 85
write(ccol,'(i4)')icol 86
write(crow,'(i4)')irow 87
C........Open the output file and write the header information. 88
open(66,file=filout,form='unformatted',status='new', 89
& access='direct',recl=1) 90
write(66,rec=1)'P' 91
write(66,rec=2)char(52+icolor) 92
write(66,rec=3)char(10) 93
write(66,rec=4)ccol(1:1) 94
write(66,rec=5)ccol(2:2) 95
write(66,rec=6)ccol(3:3) 96
write(66,rec=7)ccol(4:4) 97
write(66,rec=8)' ' 98
write(66,rec=9)crow(1:1) 99
write(66,rec=10)crow(2:2) 100
write(66,rec=11)crow(3:3) 101
write(66,rec=12)crow(4:4) 102
write(66,rec=13)char(10) 103
write(66,rec=14)'2' 104
write(66,rec=15)'5' 105
write(66,rec=16)'5' 106
write(66,rec=17)char(10) 107
C........Initialize record pointers 108
iout=18 109
iin =3 110
C........Calculate "epsilon" that is added to values if compressing. 111
if (ndec.ne.0) eps = 10.0**(-ndec)*xnorm/(1.0-10.0**(-ndec)) 112
C........Loop through all pixels in the image. 113
do 200 j=1,irow 114
do 100 i=1,icol 115
if (i .le. 10) then 116
C..............Draw the color map if within the first 10 columns. 117
if (icolor .eq. 1 .or. iside .eq. 1) then 118
xpnt = real(irow-j)/real(irow-1)*xnorm 119
else 120
xpnt = (-2./real(irow-1)*real(j) 121
& + real(irow+1)/real(irow-1))*xnorm 122
endif 123
else 124
C..............If not in first 10 columns, get data from file. 125
read(55,rec=iin,err=2000)xpnt 126
iin = iin+1 127
endif 128
C..............If the point is out of range, whack it down to size. 129
if (abs(xpnt) .gt. xnorm) xpnt = xpnt/abs(xpnt)*xnorm 130
C..............Use log scaling if requested. 131
if (ndec .ne. 0) then 132
xtmp = log10((abs(xpnt)+eps)/(xnorm+eps))/ 133
& real(ndec)+1. 134
else 135
xtmp = abs(xpnt/xnorm) 136
endif 137
if (icolor .eq. 1) then 138
C.................Grayscale output. 139
write(66,rec=iout)char(int(255*xtmp)) 140
iout = iout+1 141
else 142
if (iside .eq. 1) then 143
C.................Use color map #1 for one-sided data. 144
call cmap1(h,s,v,xtmp) 145
else 146
C.................Reconstruct sign & use color map #2 for two-sided. 147
if (xpnt .lt. 0.) xtmp = -xtmp 148
call cmap2(h,s,v,xtmp) 149
endif 150
C.................Convert (h,s,v) to (r,g,b) and write output. 151
call hsvrgb(h,s,v,ir,ig,ib) 152
write(66,rec=iout)char(ir) 153
write(66,rec=iout+1)char(ig) 154
write(66,rec=iout+2)char(ib) 155
iout = iout+3 156
endif 157
100 continue 158
200 continue 159
2000 close(55) 160
close(66) 161
goto 1000 162
end 163
C********************************************************************** 164
C.....Color mapping appropriate for one-sided data. 165
subroutine cmap1(h,s,v,xpnt) 166
real h, s, v, xpnt, pi 167
parameter (pi=3.141593) 168
if (xpnt .gt. 0.8) then 169
h = 330. + 30*(1.-xpnt)/0.2 170
s = (1.-xpnt)/0.2 171
v = 1.0 172
elseif (xpnt .gt. 0.2) then 173
t = 2.0*pi*((xpnt-0.2)/.6) 174
h = 240.0*(0.8-xpnt)/0.6 + 5.0*sin(t) 175
s = 1.0 176
v = 1.0 177
else 178
h = 270.0*(.2-xpnt)+240.0 179
s = 1.0 180
v = xpnt/0.2 181
endif 182
return 183
end 184
C********************************************************************** 185
C.....Color mapping appropriate for two-sided data. Zero = black. 186
subroutine cmap2(h,s,v,xpnt) 187
real h, s, v, xpnt 188
if (xpnt .gt. 0.8) then 189
h = 30. - 30.*(xpnt-0.8)*(xpnt-0.8)/0.04 190
s = (1.0-xpnt)/0.2 191
v = 1.0 192
else if (xpnt .gt. .3) then 193
h = 135.0 - 105.*(xpnt-.3)/.5 194
s = 1.0 195
v = 1.0 196
else if (xpnt .gt. 0.) then 197
h = 180.0 - 45.*xpnt/.3 198
s = 1.0 199
v = xpnt/.3 200
else if (xpnt .gt. -.3) then 201
h = (180. - 45.*xpnt/.3) 202
s = 1.0 203
v = abs(xpnt)/.3 204
else 205
h = 225. - 135.*(xpnt+.3)/.7 206
s = 1.0 207
v = 1.0 208
endif 209
return 210
end 211
C********************************************************************** 212
C.....Subroutine to convert (h,s,v) values to (r,g,b) values. 213
subroutine hsvrgb(h,s,v,ir,ig,ib) 214
implicit real (a-h,o-z) 215
real pqvt(4) 216
integer irpnt(6), igpnt(6), ibpnt(6) 217
data irpnt/3,2,1,1,4,3/,igpnt/4,3,3,2,1,1/,ibpnt/1,1,4,3,3,2/ 218
C..... 219
C input: 0 <= h <= 360 (hue) 220
C 0 <= v <= 1.0 (value or brightness) 221
C 0 <= s <= 1.0 (saturation) 222
C output: 0 <= ir,ig,ib <= maxval 223
C..... 224
if (s .eq. 0.0) then 225
ir = 255.0*v 226
ig = ir 227
ib = ir 228
return 229
endif 230
if (h .eq. 360.0) h=0.0 231
h = h / 60.0 232
i = int(h) 233
f = h-i 234
pqvt(3) = 255.0*v 235
pqvt(1) = pqvt(3)*(1.0-s) 236
pqvt(2) = pqvt(3)*(1.0-s*f) 237
pqvt(4) = pqvt(3)*(1.0-s*(1.0-f)) 238
i = i + 1 239
ir = pqvt(irpnt(i)) 240
ig = pqvt(igpnt(i)) 241
ib = pqvt(ibpnt(i)) 242
return 243
end 244
C********************************************************************** 245
After the input and output files names have been constructed based on the current frame number and the given base name, the status of these files is checked. If the input file is not found, the program terminates (lines 67 to 71). If the input file exists and the output file already exists, another output file is not created. Instead, the frame number is advanced and the next set of input and output files is checked (lines 73 to 77). Only if the input file exists and the output file does not exist is a new output file created. This is done in order to permit parallel processing of the data---this program may be invoked more than once to work on the same set of files and the file handling should insure that separate invocations do not interfere with each other.
The majority of code should be easy to understand by experienced Fortran programmers. However, the use of the Fortran intrinsic function char in lines 140, 153, 154, and 155 to convert an integer value (which is between 0 and 255) to the corresponding single byte value (characters are inherently one byte long) may be unfamiliar.
The subroutine cmap1 (ref. lines 165-185) maps a one-sided field to a color. The fields are assumed to be between zero and unity and the subroutine returns the corresponding location in the HSV color space.
The subroutine cmap2 (ref. lines 186-211) maps a two-sided field to a color. The fields are assumed to be between -1 and 1 and the subroutine returns the corresponding location in the HSV color space.
The subroutine hsvrgb (ref. lines 213-244) converts from HSV values to RGB values.