Add APNG backend for animations (#45)

Frames written to a .png or .apng file are stored as full-colour RGB
instead of being quantised to a 256 colour GIF palette. The deflate and
scanline filtering of png_encode is shared with the new apng_encode.

Author: jalvesz

README.md

@@ -123,7 +123,7 @@ call suptitle("figure title")
   scale along an edge
 - `table` for a block of text below (or above) the axes
 - `streamplot` for the streamlines of a vector field
-- `add_frame` / `save_animation(file, fps=, loop=)` to write an animated GIF
+- `add_frame` / `save_animation(file, fps=, loop=)` to write an animated GIF, or a full-colour animated PNG with a `.png` or `.apng` file name
 - `spy` (sparsity pattern of a dense matrix or COO indices), `matshow`, `eventplot` and `broken_barh`
 - `hist2d` and `hexbin` for counting points into square or hexagonal bins
 - `imshow(interpolation="bilinear")` to smooth an image instead of

src/fplotlib.f90

@@ -13,6 +13,7 @@ module fplotlib
     use fplotlib_backend_pdf
     use fplotlib_backend_eps
     use fplotlib_gif, only: gif_encode
+    use fplotlib_png, only: apng_encode
     use fplotlib_proj3d
     use fplotlib_backend_png
     use fplotlib_mathtext
@@ -6233,7 +6234,8 @@ subroutine anim_append(frame)
         anim_len = need
     end subroutine anim_append
 
-    ! Writes the frames collected so far and starts a new reel, so that a
+    ! Writes the frames collected so far as a GIF, or as an animated PNG for
+    ! .png and .apng (full colour), and starts a new reel, so that a
     ! program can make several animations without a reset call.
     subroutine save_animation(filename, fps, loop)
         character(len=*), intent(in) :: filename
@@ -6247,10 +6249,13 @@ subroutine save_animation(filename, fps, loop)
             print *, "fplotlib: no frames; call add_frame before save_animation"
             return
         end if
-        if (file_ext(filename) /= "gif") then
-            print *, "fplotlib: save_animation writes .gif, not .", file_ext(filename)
+        select case (file_ext(filename))
+        case ("gif", "png", "apng")
+        case default
+            print *, "fplotlib: save_animation writes .gif, .png or .apng, not .", &
+                file_ext(filename)
             return
-        end if
+        end select
         rate = 5.0_dp
         if (present(fps)) then
             if (fps <= 0.0_dp) error stop "fplotlib: save_animation fps must be positive"
@@ -6262,8 +6267,16 @@ subroutine save_animation(filename, fps, loop)
         ! GIF counts delays in hundredths of a second and viewers treat 0 as
         ! "as fast as you like", so a frame never asks for less than one.
         delay = max(1, nint(100.0_dp/rate))
-        call write_bytes(filename, gif_encode(anim_w, anim_h, anim_n, &
-                                              anim_pix(1:anim_len), delay, rep))
+        if (file_ext(filename) == "gif") then
+            call write_bytes(filename, gif_encode(anim_w, anim_h, anim_n, &
+                                                  anim_pix(1:anim_len), delay, rep))
+        else
+            ! APNG delays are a fraction, so the frame rate is not rounded
+            ! to hundredths of a second.
+            call write_bytes(filename, apng_encode(anim_w, anim_h, anim_n, &
+                             anim_pix(1:anim_len), &
+                             min(65535, max(1, nint(1000.0_dp/rate))), 1000, rep))
+        end if
         anim_n = 0
         anim_w = 0
         anim_h = 0

src/fplotlib_png.f90

@@ -21,7 +21,7 @@ module fplotlib_png
     ! zlib_compress is here rather than in the PDF backend because a PDF
     ! FlateDecode stream is the same zlib stream a PNG carries, and one
     ! compressor is enough.
-    public :: png_encode, zlib_compress
+    public :: png_encode, apng_encode, zlib_compress
 
     ! LZ77 parameters. WSIZE is fixed by the format; CHAIN bounds how hard the
     ! matcher looks before settling, which is the usual size/speed dial.
@@ -314,12 +314,113 @@ function png_encode(w, h, rgba) result(bytes)
         integer, intent(in) :: w, h
         integer, intent(in) :: rgba(4, w, h)
         character(len=:), allocatable :: bytes
+        type(bitbuf_t) :: out
+
+        call push_signature(out)
+        call ihdr(out, w, h, 6)
+        call chunk(out, "IDAT", filtered_zlib(w, h, 4, rgba))
+        call chunk(out, "IEND", NO_DATA)
+        bytes = to_string(out)
+    end function png_encode
+
+    ! ------------------------------------------------------------------
+    ! Encode nframes RGB frames as an animated PNG. pix holds the frames as
+    ! packed RGB bytes, one after another, as the animation reel keeps them.
+    ! A delay is delay_num/delay_den seconds; loop true repeats forever.
+    ! Every frame is a full, opaque replacement of the canvas, so the first
+    ! frame is also the default image that viewers without APNG support show.
+    ! ------------------------------------------------------------------
+
+    function apng_encode(w, h, nframes, pix, delay_num, delay_den, loop) result(bytes)
+        integer, intent(in) :: w, h, nframes, delay_num, delay_den
+        character(len=*), intent(in) :: pix
+        logical, intent(in) :: loop
+        character(len=:), allocatable :: bytes
+        integer, allocatable :: rgb(:, :, :), z(:), d(:)
+        integer :: f, x, y, c, k, seq, base, nz
+        type(bitbuf_t) :: out
+
+        allocate (rgb(3, w, h))
+        call push_signature(out)
+        call ihdr(out, w, h, 2)
+
+        allocate (d(8))
+        call be32(nframes, d(1:4))
+        d(5:8) = 0                  ! zero plays forever
+        if (.not. loop) call be32(1, d(5:8))
+        call chunk(out, "acTL", d)
+        deallocate (d)
+
+        seq = 0
+        do f = 1, nframes
+            base = (f - 1)*3*w*h
+            do concurrent (y = 1:h, x = 1:w, c = 1:3)
+                k = base + (y - 1)*3*w + (x - 1)*3 + c
+                rgb(c, x, y) = iachar(pix(k:k))
+            end do
+
+            allocate (d(26))
+            call be32(seq, d(1:4))
+            call be32(w, d(5:8))
+            call be32(h, d(9:12))
+            d(13:20) = 0            ! x and y offset
+            d(21) = iand(ishft(delay_num, -8), 255)
+            d(22) = iand(delay_num, 255)
+            d(23) = iand(ishft(delay_den, -8), 255)
+            d(24) = iand(delay_den, 255)
+            d(25) = 0               ! dispose: leave the canvas as it is
+            d(26) = 0               ! blend: replace
+            call chunk(out, "fcTL", d)
+            deallocate (d)
+            seq = seq + 1
+
+            z = filtered_zlib(w, h, 3, rgb)
+            if (f == 1) then
+                call chunk(out, "IDAT", z)
+            else
+                nz = size(z)
+                allocate (d(4 + nz))
+                call be32(seq, d(1:4))
+                d(5:) = z
+                call chunk(out, "fdAT", d)
+                deallocate (d)
+                seq = seq + 1
+            end if
+        end do
+
+        call chunk(out, "IEND", NO_DATA)
+        bytes = to_string(out)
+    end function apng_encode
+
+    subroutine push_signature(out)
+        type(bitbuf_t), intent(inout) :: out
+
+        call push(out, 137); call push(out, 80); call push(out, 78)
+        call push(out, 71); call push(out, 13); call push(out, 10)
+        call push(out, 26); call push(out, 10)
+    end subroutine push_signature
+
+    function to_string(out) result(bytes)
+        type(bitbuf_t), intent(in) :: out
+        character(len=:), allocatable :: bytes
+        integer :: i
+
+        allocate (character(len=out%n) :: bytes)
+        do i = 1, out%n
+            bytes(i:i) = achar(out%b(i))
+        end do
+    end function to_string
+
+    ! The zlib stream of an image's scanlines, each filtered with whichever of
+    ! the five PNG filters looks cheapest. px is bpp interleaved channels.
+    function filtered_zlib(w, h, bpp, px) result(z)
+        integer, intent(in) :: w, h, bpp
+        integer, intent(in) :: px(bpp, w, h)
+        integer, allocatable :: z(:)
         integer, allocatable :: raw(:), comp(:), cur(:), up(:), try(:), best(:)
-        integer :: bpp, stride, nraw, x, y, c, i, p, f, score, bscore, bf
+        integer :: stride, nraw, x, y, c, i, p, f, score, bscore, bf
         integer :: a_, b_, c_, v
-        type(bitbuf_t) :: out
 
-        bpp = 4
         stride = w*bpp
         nraw = h*(stride + 1)
         allocate (raw(nraw), cur(stride), up(stride), try(stride), best(stride))
@@ -329,9 +430,9 @@ function png_encode(w, h, rgba) result(bytes)
         do y = 1, h
             i = 0
             do x = 1, w
-                do c = 1, 4
+                do c = 1, bpp
                     i = i + 1
-                    cur(i) = iand(rgba(c, x, y), 255)
+                    cur(i) = iand(px(c, x, y), 255)
                 end do
             end do
 
@@ -372,21 +473,8 @@ function png_encode(w, h, rgba) result(bytes)
         end do
 
         comp = deflate(raw, nraw)
-
-        ! Signature
-        call push(out, 137); call push(out, 80); call push(out, 78)
-        call push(out, 71); call push(out, 13); call push(out, 10)
-        call push(out, 26); call push(out, 10)
-
-        call ihdr(out, w, h)
-        call chunk(out, "IDAT", zlib_wrap(comp, raw, nraw))
-        call chunk(out, "IEND", NO_DATA)
-
-        allocate (character(len=out%n) :: bytes)
-        do i = 1, out%n
-            bytes(i:i) = achar(out%b(i))
-        end do
-    end function png_encode
+        z = zlib_wrap(comp, raw, nraw)
+    end function filtered_zlib
 
     ! Bytes in, a zlib stream out. Unlike a PNG the caller supplies its own
     ! layout, so nothing is filtered here.
@@ -421,15 +509,15 @@ function zlib_wrap(comp, raw, nraw) result(z)
         z(n + 6) = iand(a, 255)
     end function zlib_wrap
 
-    subroutine ihdr(out, w, h)
+    subroutine ihdr(out, w, h, ctype)
         type(bitbuf_t), intent(inout) :: out
-        integer, intent(in) :: w, h
+        integer, intent(in) :: w, h, ctype
         integer :: d(13)
 
         call be32(w, d(1:4))
         call be32(h, d(5:8))
         d(9) = 8      ! bit depth
-        d(10) = 6     ! truecolor with alpha
+        d(10) = ctype ! 2 truecolor, 6 truecolor with alpha
         d(11) = 0     ! deflate
         d(12) = 0     ! adaptive filtering
         d(13) = 0     ! no interlace