Skip to content

Commit

Permalink
adding implcit none to jpc functions, additional testing for jpcpack/…
Browse files Browse the repository at this point in the history 
…jpcunpack (#516)

* added implicit none to jpcpack.F90

* added implicit none to jpcunpack

* more testing of jpcpack/jpcunpack

* more testing of jpcpack/jpcunpack

* more testing of jpcpack/jpcunpack

* more testing of jpcpack/jpcunpack

* more testing
  • Loading branch information
@edwardhartnett
edwardhartnett authored Aug 3, 2023
1 parent fe85cb7 commit 8e0807a
Show file tree
Hide file tree
Showing 3 changed files with 170 additions and 130 deletions.
165 changes: 81 additions & 84 deletions src/jpcpack.F90
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
!> @file
!> @brief Pack up a data field into a JPEG2000 code stream.
!> @brief Pack a data field into a JPEG2000 code stream.
!> @author Stephen Gilbert @date 2002-12-17

!> Pack up a data field into a JPEG2000 code stream.
!> Pack a data field into a JPEG2000 code stream.
!>
!> After the data field is scaled, and the reference value is
!> subtracted out, it is treated as a grayscale image and passed to a
Expand All @@ -28,13 +28,16 @@
!> the packed data in bytes.
!>
!> @author Stephen Gilbert @date 2002-12-17
subroutine jpcpack(fld,width,height,idrstmpl,cpack,lcpack)

integer,intent(in) :: width,height
real,intent(in) :: fld(width*height)
subroutine jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
implicit none

integer,intent(in) :: width, height
real,intent(in) :: fld(width * height)
character(len=1),intent(out) :: cpack(*)
integer,intent(inout) :: idrstmpl(*)
integer,intent(inout) :: lcpack
integer :: ndpts, nbytes, nbits, maxdif, j, imin, imax
real :: temp, dscale, bscale

interface
function enc_jpeg2000(cin, width, height, nbits, ltype, ratio, retry, outjpc, jpclen) &
Expand All @@ -48,111 +51,105 @@ function enc_jpeg2000(cin, width, height, nbits, ltype, ratio, retry, outjpc, jp
end function enc_jpeg2000
end interface

real(4) :: ref,rmin4
real(8) :: rmin,rmax
real(4) :: ref, rmin4
real(8) :: rmin, rmax
integer(4) :: iref
integer(8) :: nsize
integer :: ifld(width*height),retry
integer,parameter :: zero=0
character(len=1),allocatable :: ctemp(:)
integer :: ifld(width * height), retry
integer, parameter :: zero = 0
character(len = 1), allocatable :: ctemp(:)

ndpts=width*height
bscale=2.0**real(-idrstmpl(2))
dscale=10.0**real(idrstmpl(3))
ndpts = width * height
bscale = 2.0**real(-idrstmpl(2))
dscale = 10.0**real(idrstmpl(3))

! Find max and min values in the data
if(ndpts>0) then
rmax=fld(1)
rmin=fld(1)
! Find max and min values in the data.
if (ndpts > 0) then
rmax = fld(1)
rmin = fld(1)
else
rmax=1.0
rmin=1.0
rmax = 1.0
rmin = 1.0
endif
do j=2,ndpts
if (fld(j).gt.rmax) rmax=fld(j)
if (fld(j).lt.rmin) rmin=fld(j)
do j = 2, ndpts
if (fld(j) .gt. rmax) rmax = fld(j)
if (fld(j) .lt. rmin) rmin = fld(j)
enddo
if (idrstmpl(2).eq.0) then
maxdif=nint(rmax*dscale)-nint(rmin*dscale)
if (idrstmpl(2) .eq. 0) then
maxdif = nint(rmax * dscale) - nint(rmin * dscale)
else
maxdif=nint((rmax-rmin)*dscale*bscale)
maxdif = nint((rmax - rmin) * dscale * bscale)
endif

! If max and min values are not equal, pack up field.
! If they are equal, we have a constant field, and the reference
! value (rmin) is the value for each point in the field and
! set nbits to 0.
if (rmin.ne.rmax .AND. maxdif.ne.0) then

! Determine which algorithm to use based on user-supplied
! binary scale factor and number of bits.
if (idrstmpl(2).eq.0) then

! No binary scaling and calculate minimum number of
! bits in which the data will fit.
imin=nint(rmin*dscale)
imax=nint(rmax*dscale)
maxdif=imax-imin
temp=alog(real(maxdif+1))/alog(2.0)
nbits=ceiling(temp)
rmin=real(imin)
! scale data
do j=1,ndpts
ifld(j)=nint(fld(j)*dscale)-imin
! If max and min values are not equal, pack up field. If they are
! equal, we have a constant field, and the reference value (rmin) is
! the value for each point in the field and set nbits to 0.
if (rmin .ne. rmax .AND. maxdif .ne. 0) then
! Determine which algorithm to use based on user-supplied binary
! scale factor and number of bits.
if (idrstmpl(2) .eq. 0) then
! No binary scaling and calculate minimum number of bits in
! which the data will fit.
imin = nint(rmin * dscale)
imax = nint(rmax * dscale)
maxdif = imax - imin
temp = alog(real(maxdif + 1)) / alog(2.0)
nbits = ceiling(temp)
rmin = real(imin)
! Scale data.
do j = 1, ndpts
ifld(j) = nint(fld(j) * dscale) - imin
enddo
else

! Use binary scaling factor and calculate minimum number of
! bits in which the data will fit.
rmin=rmin*dscale
rmax=rmax*dscale
maxdif=nint((rmax-rmin)*bscale)
temp=alog(real(maxdif+1))/alog(2.0)
nbits=ceiling(temp)
! Use binary scaling factor and calculate minimum number of
! bits in which the data will fit.
rmin = rmin * dscale
rmax = rmax * dscale
maxdif = nint((rmax - rmin) * bscale)
temp = alog(real(maxdif + 1)) / alog(2.0)
nbits = ceiling(temp)
! scale data
do j=1,ndpts
ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
do j = 1, ndpts
ifld(j) = max(0, nint(((fld(j) * dscale) - rmin) * bscale))
enddo
endif

! Pack data into full octets, then do JPEG2000 encode.
! and calculate the length of the packed data in bytes
retry=0
nbytes=(nbits+7)/8
nsize=lcpack ! needed for input to enc_jpeg2000
allocate(ctemp(nbytes*ndpts))
call g2_sbytesc(ctemp,ifld,0,nbytes*8,0,ndpts)
lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6), &
idrstmpl(7),retry,cpack,nsize)
if (lcpack.le.0) then
print *,'jpcpack: ERROR Packing JPC=',lcpack
if (lcpack.eq.-3) then
retry=1
! Pack data into full octets, then do JPEG2000 encode. and
! calculate the length of the packed data in bytes
retry = 0
nbytes = (nbits + 7) / 8
nsize = lcpack ! needed for input to enc_jpeg2000
allocate(ctemp(nbytes * ndpts))
call g2_sbytesc(ctemp, ifld, 0, nbytes * 8, 0, ndpts)
lcpack = enc_jpeg2000(ctemp, width, height, nbits, idrstmpl(6), &
idrstmpl(7), retry, cpack, nsize)
if (lcpack .le. 0) then
print *,'jpcpack: ERROR Packing JPC = ',lcpack
if (lcpack .eq. -3) then
retry = 1
print *,'jpcpack: Retrying....'
lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6), &
idrstmpl(7),retry,cpack,nsize)
if (lcpack.le.0) then
lcpack = enc_jpeg2000(ctemp, width, height, nbits, idrstmpl(6), &
idrstmpl(7), retry, cpack, nsize)
if (lcpack .le. 0) then
print *,'jpcpack: Retry Failed.'
else
print *,'jpcpack: Retry Successful.'
endif
endif
endif
deallocate(ctemp)

else
nbits=0
lcpack=0
nbits = 0
lcpack = 0
endif

! Fill in ref value and number of bits in Template 5.0
! Fill in ref value and number of bits in Template 5.0.
rmin4 = real(rmin, 4)
call mkieee(rmin4,ref,1) ! ensure reference value is IEEE format
! call g2_gbytec(ref,idrstmpl(1),0,32)
iref=transfer(ref,iref)
idrstmpl(1)=iref
idrstmpl(4)=nbits
idrstmpl(5)=0 ! original data were reals
if (idrstmpl(6).eq.0) idrstmpl(7)=255 ! lossy not used
call mkieee(rmin4, ref, 1) ! ensure reference value is IEEE format.
iref = transfer(ref, iref)
idrstmpl(1) = iref
idrstmpl(4) = nbits
idrstmpl(5) = 0 ! original data were reals.
if (idrstmpl(6) .eq. 0) idrstmpl(7) = 255 ! lossy not used

end subroutine
82 changes: 42 additions & 40 deletions src/jpcunpack.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,48 +16,50 @@
!> @param[out] fld Contains the unpacked data values.
!>
!> @author Stephen Gilbert @date 2002-12-17
subroutine jpcunpack(cpack,len,idrstmpl,ndpts,fld)
subroutine jpcunpack(cpack,len,idrstmpl,ndpts,fld)
implicit none

character(len=1),intent(in) :: cpack(len)
integer,intent(in) :: ndpts,len
integer,intent(in) :: idrstmpl(*)
real,intent(out) :: fld(ndpts)

character(len=1),intent(in) :: cpack(len)
integer,intent(in) :: ndpts,len
integer,intent(in) :: idrstmpl(*)
real,intent(out) :: fld(ndpts)
integer :: ifld(ndpts)
integer(4) :: ieee
integer(8) :: len8
real :: ref,bscale,dscale
integer :: nbits, j, iret

integer :: ifld(ndpts)
integer(4) :: ieee
integer(8) :: len8
real :: ref,bscale,dscale
interface
function dec_jpeg2000(cin, len, ifld) &
bind(c, name="g2c_dec_jpeg2000")
use iso_c_binding
character(kind = c_char), intent(in) :: cin(*)
integer(c_size_t), value, intent(in) :: len
integer(c_int), intent(inout) :: ifld(*)
integer(c_int) :: dec_jpeg2000
end function dec_jpeg2000
end interface

interface
function dec_jpeg2000(cin, len, ifld) &
bind(c, name="g2c_dec_jpeg2000")
use iso_c_binding
character(kind = c_char), intent(in) :: cin(*)
integer(c_size_t), value, intent(in) :: len
integer(c_int), intent(inout) :: ifld(*)
integer(c_int) :: dec_jpeg2000
end function dec_jpeg2000
end interface
ieee = idrstmpl(1)
call rdieee(ieee,ref,1)
bscale = 2.0**real(idrstmpl(2))
dscale = 10.0**real(-idrstmpl(3))
nbits = idrstmpl(4)

ieee = idrstmpl(1)
call rdieee(ieee,ref,1)
bscale = 2.0**real(idrstmpl(2))
dscale = 10.0**real(-idrstmpl(3))
nbits = idrstmpl(4)
! if nbits equals 0, we have a constant field where the reference value
! is the data value at each gridpoint
if (nbits.ne.0) then
! call g2_gbytesc(cpack,ifld,0,nbits,0,ndpts)
len8 = len
iret=dec_jpeg2000(cpack,len8,ifld)
do j=1,ndpts
fld(j)=((real(ifld(j))*bscale)+ref)*dscale
enddo
else
do j=1,ndpts
fld(j)=ref
enddo
endif

! if nbits equals 0, we have a constant field where the reference value
! is the data value at each gridpoint
if (nbits.ne.0) then
! call g2_gbytesc(cpack,ifld,0,nbits,0,ndpts)
len8 = len
iret=dec_jpeg2000(cpack,len8,ifld)
do j=1,ndpts
fld(j)=((real(ifld(j))*bscale)+ref)*dscale
enddo
else
do j=1,ndpts
fld(j)=ref
enddo
endif

end
end subroutine jpcunpack
53 changes: 47 additions & 6 deletions tests/test_jpcpack.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
! project.
!
! Brian Curtis 2021-10-18
! Ed Hartnett, 7/31/23
program test_jpcpack

implicit none

integer, parameter :: width = 2, height = 2, ndpts = 4
Expand All @@ -27,23 +27,64 @@ program test_jpcpack
idrstmpl(6) = 0
idrstmpl(7) = 1

print *, 'Testing jpcpack/jpcunpack with no data...'

! Pack the data.
call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
call jpcpack(fld, 0, height, idrstmpl, cpack, lcpack)
print *, 'lcpack: ', lcpack
if (lcpack .ne. 0) stop 1

! Unpack the data.
call jpcunpack(cpack, lcpack, idrstmpl, 0, fld2)

! Compare each value to see match, remember, comparing reals.
print *, 'ok!'
print *, 'Testing jpcpack/jpcunpack with data...'

! Pack the data.
lcpack = 200
call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
! The size of lcpack will depend on the version of jasper used to compress the data.
!print *,lcpack
if (lcpack .le. 0) stop 2

! Unpack the data.
call jpcunpack(cpack, lcpack, idrstmpl, ndpts, fld2)

! Compare each value to see match, remember, comparing reals
print *, fld
print *, fld2
! Compare each value to see match, remember, comparing reals.
!print *, fld
!print *, fld2
do ii = 1, ndpts
if (abs(fld(ii) - fld2(ii)) .gt. delta) then
print *, fld(ii), fld2(ii), 'do not match'
stop 4
end if
end do

print *, 'SUCCESS!'
print *, 'ok!'
print *, 'Testing jpcpack/jpcunpack with data and idrstmpl(2) = 0...'

! Pack the data.
idrstmpl(2) = 0
lcpack = 200
call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
! The size of lcpack will depend on the version of jasper used to compress the data.
!print *,lcpack
if (lcpack .le. 0) stop 2

! Unpack the data.
call jpcunpack(cpack, lcpack, idrstmpl, ndpts, fld2)

! Compare each value to see match, remember, comparing reals.
!print *, fld
!print *, fld2
do ii = 1, ndpts
if (abs(fld(ii) - fld2(ii)) .gt. delta) then
print *, fld(ii), fld2(ii), 'do not match'
stop 4
end if
end do

print *, 'ok!'
print *, 'SUCCESS!'
end program test_jpcpack

0 comments on commit 8e0807a

Please sign in to comment.