;;;; Apress-fractal.lisp -- ;;;; An entry into the http://www.apress.com/promo/fractal/ contest ;;;; ;;;; Copyright (c) 2005 ;;;; ;;;; By David Steuber ;;;; ;;;; Permission is hereby granted, free of charge, to any person obtaining ;;;; a copy of this software and associated documentation files (the ;;;; "Software"), to deal in the Software without restriction, including ;;;; without limitation the rights to use, copy, modify, merge, publish, ;;;; distribute, sublicense, and/or sell copies of the Software, and to ;;;; permit persons to whom the Software is furnished to do so, subject to ;;;; the following conditions: ;;;; ;;;; The above copyright notice and this permission notice shall be ;;;; included in all copies or substantial portions of the Software. ;;;; ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF ;;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE ;;;; LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION ;;;; OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION ;;;; WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ;;;; ;;;; *********************************************************************** ;;;; *********************************************************************** ;;;; *********************************************************************** ;;;; *********************************************************************** ;;;; *********************************************************************** ;;;; ;;;; JPEG Encoder (edited by David Steuber to reduce functionality ;;;; and those pesky, but useful, optimization notes from SBCL): ;;;; ;;;; Copyright (c) 1999, Eugene Zaikonnikov ;;;; All rights reserved. ;;;; ;;;; Redistribution and use in source and binary forms, with or without ;;;; modification, are permitted provided that the following conditions are ;;;; met: ;;;; ;;;; * Redistributions of source code must retain the above copyright ;;;; notice, this list of conditions and the following disclaimer. ;;;; ;;;; * Redistributions in binary form must reproduce the above ;;;; copyright notice, this list of conditions and the following disclaimer ;;;; in the documentation and/or other materials provided with the ;;;; distribution. ;;;; ;;;; * Neither the name of the author nor the names of its contributors ;;;; may be used to endorse or promote products derived from this software ;;;; without specific prior written permission. ;;;; ;;;; THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ;;;; "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ;;;; LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ;;;; A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ;;;; OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ;;;; SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ;;;; LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, ;;;; DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ;;;; THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;;;; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ;;;; OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;;;; ************************************************************************ ;;; Generic Common Lisp JPEG encoder/decoder implementation ;;; $Id: jpeg.lisp,v 1.8 2003/10/06 17:49:30 kevinrosenberg Exp $ ;;; Version 1.022, June 1999. ;;; Written by Eugene Zaikonnikov [viking@funcall.org] ;;; Copyright [c] 1999, Eugene Zaikonnikov ;;; This software is distributed under the terms of BSD-like license ;;; [see LICENSE for details] ;;; That was qute some time ago - I'd wrote it better now [E.Z., 2001] ;;; Known to work with Lispworks 4 and Allegro CL 5 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Creation of this software was sponsored by Kelly E. Murray ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; WARNING: IT IS HIGHLY RECOMMENDED TO SUSPEND ALL CPU-CONSUMING OS THREADS DURING COMPILATION. OTHERWISE SYSTEM-DEPENDENT ;;; OPTIMISATIONS MAY BE MISIDENTIFIED WHICH SHOULD RESULT IN A SUFFICIENTLY DEGRADED PERFORMANCE. ;;; Two main functions available: ;;; ;;; (encode-image filename image ncomp h w &key sampling q-tabs q-factor), where: ;;; filename - output file name ;;; ncomp - number of components (1-4) ;;; h, w - source image height and width respectively ;;; image - array of B, G, R pixels in case of three component image, array of grayscale pixels in case of single component, ;;; array of 2 or 4 pixels in the case of two or four component image respectively ;;; :q-tabs - specifies quantization tables vector, should be 1 for 1, 2 for 2, 2 for 3 and 4 entries for 4 components ;;; :sampling - sampling frequency for ncomp components by X and Y axis, e.g. '((2 2) (1 1) (1 1)) for three components, can be omitted ;;; for grayscale and RGB images ;;; :q-factor - quality specifier (1-64), default is 64 ;;; Returns nothing of practical use ;;; ;;; Based on CCITT Rec. T.81 "Information technology - digital compression and coding of continious-tone still images - ;;; requirements and guidelines". ;;; Credits: ;;; to the Independent JPEG Group - colorspace conversion and DCT algorithms were adopted from their sources; ;;; to Jeff Dalton for his wise paper "Common Lisp Pitfalls". (defpackage #:jpeg (:use #:common-lisp)) (in-package #:jpeg) (eval-when (:compile-toplevel) (export '(encode-image))) (declaim (inline csize write-stuffed quantize get-average zigzag encode-block llm-dct descale crunch colorspace-convert subsample inverse-llm-dct dequantize upsample extend recieve decode-ac decode-dc decode-block izigzag write-bits decompose-yuv)) (eval-when (:compile-toplevel :load-toplevel :execute) (defvar *optimize* '(optimize (safety 0) (space 0) (debug 0) (speed 3)))) ; '(optimize (safety 1) (space 1) (debug 1) (speed 0)))) (eval-when (:compile-toplevel :load-toplevel :execute) ;;; For ease of reference (defmacro dbref (data x y) `(the fixnum (svref (svref ,data ,y) ,x))) ;;; Integer arithmetic wrappers (defmacro plus (a b) `(the fixnum (+ (the fixnum ,a) (the fixnum ,b)))) (defmacro minus (a b) `(the fixnum (- (the fixnum ,a) (the fixnum ,b)))) (defmacro mul (a b) `(the fixnum (* (the fixnum ,a) (the fixnum ,b)))) ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Here we define some constants (markers, quantization and huffman tables etc.) (eval-when (:compile-toplevel :load-toplevel) (defmacro defconstant* (name value &optional doc) "Ensure VALUE is evaluated only once, required for ANSI and SBCL." `(defconstant ,name (if (boundp ',name) (symbol-value ',name) ,value) ,@(when doc (list doc)))) ;;; Most positive (signed-byte 24) (defconstant* +most-positive-signed-byte-24+ (1- (ash 1 23))) ;;; Source huffman tables for the encoder (defconstant* +luminance-dc-bits+ #(#x00 #x01 #x05 #x01 #x01 #x01 #x01 #x01 #x01 #x00 #x00 #x00 #x00 #x00 #x00 #x00)) (defconstant* +luminance-dc-values+ #(#x00 #x01 #x02 #x03 #x04 #x05 #x06 #x07 #x08 #x09 #x0a #x0b)) (defconstant* +chrominance-dc-bits+ #(#x00 #x03 #x01 #x01 #x01 #x01 #x01 #x01 #x01 #x01 #x01 #x00 #x00 #x00 #x00 #x00)) (defconstant* +chrominance-dc-values+ #(#x00 #x01 #x02 #x03 #x04 #x05 #x06 #x07 #x08 #x09 #x0a #x0b)) (defconstant* +luminance-ac-bits+ #(#x00 #x02 #x01 #x03 #x03 #x02 #x04 #x03 #x05 #x05 #x04 #x04 #x00 #x00 #x01 #x7d)) (defconstant* +luminance-ac-values+ #(#x01 #x02 #x03 #x00 #x04 #x11 #x05 #x12 #x21 #x31 #x41 #x06 #x13 #x51 #x61 #x07 #x22 #x71 #x14 #x32 #x81 #x91 #xa1 #x08 #x23 #x42 #xb1 #xc1 #x15 #x52 #xd1 #xf0 #x24 #x33 #x62 #x72 #x82 #x09 #x0a #x16 #x17 #x18 #x19 #x1a #x25 #x26 #x27 #x28 #x29 #x2a #x34 #x35 #x36 #x37 #x38 #x39 #x3a #x43 #x44 #x45 #x46 #x47 #x48 #x49 #x4a #x53 #x54 #x55 #x56 #x57 #x58 #x59 #x5a #x63 #x64 #x65 #x66 #x67 #x68 #x69 #x6a #x73 #x74 #x75 #x76 #x77 #x78 #x79 #x7a #x83 #x84 #x85 #x86 #x87 #x88 #x89 #x8a #x92 #x93 #x94 #x95 #x96 #x97 #x98 #x99 #x9a #xa2 #xa3 #xa4 #xa5 #xa6 #xa7 #xa8 #xa9 #xaa #xb2 #xb3 #xb4 #xb5 #xb6 #xb7 #xb8 #xb9 #xba #xc2 #xc3 #xc4 #xc5 #xc6 #xc7 #xc8 #xc9 #xca #xd2 #xd3 #xd4 #xd5 #xd6 #xd7 #xd8 #xd9 #xda #xe1 #xe2 #xe3 #xe4 #xe5 #xe6 #xe7 #xe8 #xe9 #xea #xf1 #xf2 #xf3 #xf4 #xf5 #xf6 #xf7 #xf8 #xf9 #xfa)) (defconstant* +chrominance-ac-bits+ #(#x00 #x02 #x01 #x02 #x04 #x04 #x03 #x04 #x07 #x05 #x04 #x04 #x00 #x01 #x02 #x77)) (defconstant* +chrominance-ac-values+ #(#x00 #x01 #x02 #x03 #x11 #x04 #x05 #x21 #x31 #x06 #x12 #x41 #x51 #x07 #x61 #x71 #x13 #x22 #x32 #x81 #x08 #x14 #x42 #x91 #xa1 #xb1 #xc1 #x09 #x23 #x33 #x52 #xf0 #x15 #x62 #x72 #xd1 #x0a #x16 #x24 #x34 #xe1 #x25 #xf1 #x17 #x18 #x19 #x1a #x26 #x27 #x28 #x29 #x2a #x35 #x36 #x37 #x38 #x39 #x3a #x43 #x44 #x45 #x46 #x47 #x48 #x49 #x4a #x53 #x54 #x55 #x56 #x57 #x58 #x59 #x5a #x63 #x64 #x65 #x66 #x67 #x68 #x69 #x6a #x73 #x74 #x75 #x76 #x77 #x78 #x79 #x7a #x82 #x83 #x84 #x85 #x86 #x87 #x88 #x89 #x8a #x92 #x93 #x94 #x95 #x96 #x97 #x98 #x99 #x9a #xa2 #xa3 #xa4 #xa5 #xa6 #xa7 #xa8 #xa9 #xaa #xb2 #xb3 #xb4 #xb5 #xb6 #xb7 #xb8 #xb9 #xba #xc2 #xc3 #xc4 #xc5 #xc6 #xc7 #xc8 #xc9 #xca #xd2 #xd3 #xd4 #xd5 #xd6 #xd7 #xd8 #xd9 #xda #xe2 #xe3 #xe4 #xe5 #xe6 #xe7 #xe8 #xe9 #xea #xf2 #xf3 #xf4 #xf5 #xf6 #xf7 #xf8 #xf9 #xfa)) ;;;Zigzag encoding matrix (defconstant* +zigzag-index+ #(#(0 1 5 6 14 15 27 28) #(2 4 7 13 16 26 29 42) #(3 8 12 17 25 30 41 43) #(9 11 18 24 31 40 44 53) #(10 19 23 32 39 45 52 54) #(20 22 33 38 46 51 55 60) #(21 34 37 47 50 56 59 61) #(35 36 48 49 57 58 62 63))) ;;; Temporary buffer for zigzag encoding and decoding (defvar *zz-result* (make-array 64 :element-type 'unsigned-byte)) ;;;JPEG file markers (defconstant +M_COM+ #xfe) (defconstant +M_SOF0+ #xc0) (defconstant +M_DHT+ #xc4) (defconstant +M_RST0+ #xd0) (defconstant +M_RST7+ #xd7) (defconstant +M_SOI+ #xd8) (defconstant +M_EOI+ #xd9) (defconstant +M_SOS+ #xda) (defconstant +M_DQT+ #xdb) (defconstant +M_DNL+ #xdc) (defconstant +M_DRI+ #xdd) (defconstant +M_DAC+ #xcc) (defconstant +M_APP0+ #xe0) ;;; Default quantization tables (defvar *q-luminance* #(#(16 11 10 16 24 40 51 61) #(12 12 14 19 26 58 60 55) #(14 13 16 24 40 57 69 56) #(14 17 22 29 51 87 80 62) #(18 22 37 56 68 109 103 77) #(24 35 55 64 81 104 113 92) #(49 64 78 87 103 121 120 101) #(72 92 95 98 112 100 103 99))) (defvar *q-chrominance* #(#(17 18 24 47 99 99 99 99) #(18 21 26 66 99 99 99 99) #(24 26 56 99 99 99 99 99) #(47 66 99 99 99 99 99 99) #(99 99 99 99 99 99 99 99) #(99 99 99 99 99 99 99 99) #(99 99 99 99 99 99 99 99) #(99 99 99 99 99 99 99 99))) (defvar *q-luminance-hi* #(#(10 7 6 10 15 25 32 38) #(8 8 9 12 16 36 38 34) #(9 8 10 15 25 36 43 35) #(9 11 14 18 32 54 50 39) #(11 14 23 35 42 68 64 48) #(15 22 34 40 51 65 71 58) #(31 40 49 54 64 76 75 63) #(45 58 59 61 70 62 64 62))) (defvar *q-chrominance-hi* #(#(11 11 15 29 62 62 62 62) #(11 13 16 41 62 62 62 62) #(15 16 35 62 62 62 62 62) #(29 41 62 62 62 62 62 62) #(62 62 62 62 62 62 62 62) #(62 62 62 62 62 62 62 62) #(62 62 62 62 62 62 62 62) #(62 62 62 62 62 62 62 62))) ) ;;; Quantization performance test, each branch quantizes 3000 random matrixes ;;; Note: if your system performs this test faster than get-internal-run-time quantum, then it doesn't matters which to use (eval-when (:compile-toplevel) (format t "Performing compile-time optimization.. please wait.~%") (finish-output) (defvar *quantize-optimization* (<= (let ((time1 (get-internal-run-time))) (loop for i fixnum from 1 to 3000 do (loop for row across *q-luminance* do (loop for q-coef fixnum across row maximize (round (random 128) q-coef)))) (- (get-internal-run-time) time1)) (let ((time1 (get-internal-run-time))) (loop for i fixnum from 1 to 3000 do (loop for q-row across *q-luminance* do (loop for val fixnum = (random 128) for absval fixnum = (abs val) for qc fixnum across q-row maximize (cond ((< absval (ash qc -1)) 0) ((<= absval qc) (if (minusp val) -1 1)) ((<= (ash absval -1) qc) (if (zerop (logand absval 1)) (if (minusp val) -1 1) (if (minusp val) -2 2))) (t (round val qc)))))) (- (get-internal-run-time) time1)))) (format t "Done.~%") (finish-output) ) (defconstant* +q-tables+ (vector *q-luminance* *q-chrominance*)) ;;; This table is used to map coefficients into SSSS value (defconstant* +csize+ (make-array 2047 :initial-contents (loop for i fixnum from 0 to 2046 collecting (integer-length (abs (minus i 1023)))))) ;;; Some constants for colorspace mapper (defconstant +shift+ (1- (integer-length (ash most-positive-fixnum -7)))) (defconstant +.299+ (round (+ (* 0.299 (ash 1 +shift+)) 0.5))) (defconstant +.587+ (round (+ (* 0.587 (ash 1 +shift+)) 0.5))) (defconstant +.114+ (round (+ (* 0.114 (ash 1 +shift+)) 0.5))) (defconstant +-.1687+ (round (+ (* -0.1687 (ash 1 +shift+)) 0.5))) (defconstant +-.3313+ (round (+ (* -0.3313 (ash 1 +shift+)) 0.5))) (defconstant +-.4187+ (round (+ (* -0.4187 (ash 1 +shift+)) 0.5))) (defconstant +-.0813+ (round (+ (* -0.0813 (ash 1 +shift+)) 0.5))) (defconstant +.5+ (round (+ (* 0.5 (ash 1 +shift+)) 0.5))) (defconstant +uvoffset+ (ash 128 +shift+)) (defconstant +onehalf+ (1- (ash 1 (1- +shift+)))) (defconstant +r-y-off+ 0) (defconstant +g-y-off+ 256) (defconstant +b-y-off+ (* 2 256)) (defconstant +r-u-off+ (* 3 256)) (defconstant +g-u-off+ (* 4 256)) (defconstant +b-u-off+ (* 5 256)) (defconstant +r-v-off+ +b-u-off+) (defconstant +g-v-off+ (* 6 256)) (defconstant +b-v-off+ (* 7 256)) ;;;Direct color conversion table (defvar *ctab* (make-array 2048 :initial-element 0)) ;;; Filling in the table (loop for i fixnum from 0 to 255 do (setf (svref *ctab* (plus i +r-y-off+)) (mul +.299+ i)) (setf (svref *ctab* (plus i +g-y-off+)) (mul +.587+ i)) (setf (svref *ctab* (plus i +b-y-off+)) (mul +.114+ i)) (setf (svref *ctab* (plus i +r-u-off+)) (mul +-.1687+ i)) (setf (svref *ctab* (plus i +g-u-off+)) (mul +-.3313+ i)) (setf (svref *ctab* (plus i +b-u-off+)) (+ (mul +.5+ i) +uvoffset+ +onehalf+)) (setf (svref *ctab* (plus i +r-v-off+)) (+ (mul +.5+ i) +uvoffset+ +onehalf+)) (setf (svref *ctab* (plus i +g-v-off+)) (mul +-.4187+ i)) (setf (svref *ctab* (plus i +b-v-off+)) (mul +-.0813+ i))) ;;; Constantsants for the inverse colorspace conversion (defconstant +1.40200+ (round (+ (* 1.40200 (ash 1 +shift+)) 0.5))) (defconstant +1.77200+ (round (+ (* 1.77200 (ash 1 +shift+)) 0.5))) (defconstant +-0.71414+ (round (+ (* -0.71414 (ash 1 +shift+)) 0.5))) (defconstant +-0.34414+ (round (+ (* -0.34414 (ash 1 +shift+)) 0.5))) ;;; Inverse color conversion tables (defvar *cr-r-tab* (make-array 256)) (defvar *cb-g-tab* (make-array 256)) (defvar *cr-g-tab* (make-array 256)) (defvar *cb-b-tab* (make-array 256)) ;;; Filling up the tables (loop for i from 0 to 255 for x from -127 do (setf (svref *cr-r-tab* i) (ash (plus (mul +1.40200+ x) +onehalf+) (- +shift+))) (setf (svref *cb-b-tab* i) (ash (plus (mul +1.77200+ x) +onehalf+) (- +shift+))) (setf (svref *cr-g-tab* i) (mul +-0.71414+ x)) (setf (svref *cb-g-tab* i) (plus (mul +-0.34414+ x) +onehalf+))) ;;; Temporary workspace for IDCT (defvar *ws* (make-array 8 :initial-contents (loop for i from 0 to 7 collecting (make-array 8)))) ;;; Constants for LLM DCT (defconstant +dct-shift+ 11) ;defining DCT scaling (defconstant +shift-1+ (1- +dct-shift+)) (defconstant +shift+1+ (1+ +dct-shift+)) (defconstant +shift+4+ (+ +dct-shift+ 4)) (defconstant +FIX-0-298631336+ (round (+ (* 0.298631336 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-0-390180644+ (round (+ (* 0.390180644 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-0-541196100+ (round (+ (* 0.541196100 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-0-765366865+ (round (+ (* 0.765366865 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-0-899976223+ (round (+ (* 0.899976223 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-1-175875602+ (round (+ (* 1.175875602 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-1-501321110+ (round (+ (* 1.501321110 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-1-847759065+ (round (+ (* 1.847759065 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-1-961570560+ (round (+ (* 1.961570560 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-2-053119869+ (round (+ (* 2.053119869 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-2-562915447+ (round (+ (* 2.562915447 (ash 1 +dct-shift+)) 0.5))) (defconstant +FIX-3-072711026+ (round (+ (* 3.072711026 (ash 1 +dct-shift+)) 0.5))) ;;; Post-IDCT limiting array (defvar *idct-limit-array* (make-array 512 :initial-element 0)) (loop for n from 0 for i from 128 to 383 do (setf (svref *idct-limit-array* i) n)) (loop for i from 384 to 511 do (setf (svref *idct-limit-array* i) 255)) ;;; State variables for write-bits (defvar *prev-byte* 0) (defvar *prev-length* 0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Encoder part ;;; Subsamples inbuf into outbuf (defun subsample (inbuf outbuf H V xlim ylim iH iV) (declare #.*optimize* (type fixnum H V xlim ylim iV iH) (type (simple-vector *) inbuf outbuf)) (loop for by fixnum from 0 below V do (loop for bx fixnum from 0 below H for block = (svref outbuf (plus bx (mul by H))) do (loop for y fixnum from (ash by 3) by iV for yp fixnum from 0 to 7 do (loop for x fixnum from (ash bx 3) by iH for xp fixnum from 0 to 7 do (setf (dbref block xp yp) (the fixnum (cond ((and (<= x xlim) (<= y ylim)) (dbref inbuf x y)) ((and (> x xlim) (> y ylim)) (dbref inbuf xlim ylim)) ((> x xlim) (dbref inbuf xlim y)) ((> y ylim) (dbref inbuf x ylim)) (t (error "Internal error")) )))))))) ;;; Returns sum of Vi*Hi (defun count-relation (smp) (loop for entry in smp summing (mul (first entry) (second entry)))) ;;; Direct color mapping (defun decompose-yuv (comp-yuv YUV dx dy h w height width y-comp-lut u-comp-lut v-comp-lut) (let ((xend (plus dx (1- width))) (yend (plus dy (1- height))) (Y (aref YUV 0)) (U (aref YUV 1)) (V (aref YUV 2))) (declare #.*optimize* (type fixnum dx dy h w height width xend yend) (type (simple-vector *) y-comp-lut u-comp-lut v-comp-lut YUV Y U V) (type (simple-array (unsigned-byte 8) *) comp-yuv)) (setf xend (min xend (1- w))) (setf yend (min yend (1- h))) (loop for yd fixnum from dy to yend for ypos fixnum = (* w yd) do (loop for xd fixnum from dx to xend for pos fixnum = (plus xd ypos) for comp = (aref comp-yuv pos) for v-comp fixnum = (svref v-comp-lut comp) for u-comp fixnum = (svref u-comp-lut comp) for y-comp fixnum = (svref y-comp-lut comp) for cx fixnum = (minus xd dx) for cy fixnum = (minus yd dy) do (setf (dbref Y cx cy) y-comp) (setf (dbref U cx cy) u-comp) (setf (dbref V cx cy) v-comp))) (values xend yend))) (defun y-comp (r g b) (minus (ash (+ (svref *ctab* (plus r +r-y-off+)) (svref *ctab* (plus g +g-y-off+)) (svref *ctab* (plus b +b-y-off+))) (- +shift+)) 128)) (defun u-comp (r g b) (minus (ash (+ (svref *ctab* (plus r +r-u-off+)) (svref *ctab* (plus g +g-u-off+)) (svref *ctab* (plus b +b-u-off+))) (- +shift+)) 128)) (defun v-comp (r g b) (minus (ash (+ (svref *ctab* (plus r +r-v-off+)) (svref *ctab* (plus g +g-v-off+)) (svref *ctab* (plus b +b-v-off+))) (- +shift+)) 128)) ;;; Converts given image sampling into frequencies of pixels of components (defun convert-sampling (s Hmax Vmax) (declare (type fixnum Hmax Vmax)) (make-array (length s) :initial-contents (loop for entry in s collecting (list (the fixnum (/ Hmax (first entry))) (the fixnum (/ Vmax (second entry))))))) ;;; Quantization (also removes factor of 8 after DCT) (defmacro quantize-block () (if *quantize-optimization* '(loop for block-row across block for q-row across q-table do (loop for x fixnum from 0 to 7 for val fixnum = (ash (the fixnum (svref block-row x)) -3) for qc fixnum = (svref q-row x) do (setf (svref block-row x) (the fixnum (round val qc))))) '(loop for block-row across block for q-row across q-table do (loop for x fixnum from 0 to 7 for val fixnum = (ash (the fixnum (svref block-row x)) -3) for absval fixnum = (abs val) for qc fixnum = (svref q-row x) do (cond ((< absval (ash qc -1)) ;you won't believe, but under LWW 4.1 such ugly hack gives (setf (svref block-row x) 0)) ;very sufficient speedup ((<= absval qc) (if (minusp val) (setf (svref block-row x) -1) (setf (svref block-row x) 1))) ((<= (ash absval -1) qc) (if (zerop (logand absval 1)) (if (minusp val) (setf (svref block-row x) -1) (setf (svref block-row x) 1)) (if (minusp val) (setf (svref block-row x) -2) (setf (svref block-row x) 2)))) (t (setf (svref block-row x) (the fixnum (round val qc))))))))) (defun quantize (block q-table) (declare #.*optimize* (type (simple-vector *) block q-table)) (quantize-block)) ;;; LLM DCT aux definitions (defun descale (x n) (declare #.*optimize* (type fixnum x n)) (the fixnum (ash (plus x (ash 1 (1- n))) (- n)))) (defmacro plus3 (x y z) `(plus (plus ,x ,y) ,z)) ;;; Implementation of Loeffer, Ligtenberg and Moschytz forward DCT (defun llm-dct (data) (declare #.*optimize* (type (simple-vector *) data)) (loop with tmp0 fixnum and tmp1 fixnum and tmp2 fixnum and tmp3 fixnum and tmp4 fixnum and tmp5 fixnum and tmp6 fixnum and tmp7 fixnum and tmp10 fixnum and tmp11 fixnum and tmp12 fixnum and tmp13 fixnum and z1 fixnum and z2 fixnum and z3 fixnum and z4 fixnum and z5 fixnum do (loop ;for dptrpos fixnum from 7 downto 0 ;for dptr = (svref data dptrpos) do for dptr across data do ;iterating over rows (setf tmp0 (plus (svref dptr 0) (svref dptr 7))) (setf tmp7 (minus (svref dptr 0) (svref dptr 7))) (setf tmp1 (plus (svref dptr 1) (svref dptr 6))) (setf tmp6 (minus (svref dptr 1) (svref dptr 6))) (setf tmp2 (plus (svref dptr 2) (svref dptr 5))) (setf tmp5 (minus (svref dptr 2) (svref dptr 5))) (setf tmp3 (plus (svref dptr 3) (svref dptr 4))) (setf tmp4 (minus (svref dptr 3) (svref dptr 4))) (setf tmp10 (plus tmp0 tmp3)) (setf tmp13 (minus tmp0 tmp3)) (setf tmp11 (plus tmp1 tmp2)) (setf tmp12 (minus tmp1 tmp2)) (setf (svref dptr 0) (ash (plus tmp10 tmp11) 1)) (setf (svref dptr 4) (ash (minus tmp10 tmp11) 1)) (setf z1 (mul (plus tmp12 tmp13) +fix-0-541196100+)) (setf (svref dptr 2) (descale (plus z1 (mul tmp13 +fix-0-765366865+)) +shift-1+)) (setf (svref dptr 6) (descale (plus z1 (mul tmp12 (- +fix-1-847759065+))) +shift-1+)) (setf z1 (plus tmp4 tmp7)) (setf z2 (plus tmp5 tmp6)) (setf z3 (plus tmp4 tmp6)) (setf z4 (plus tmp5 tmp7)) (setf z5 (mul (plus z3 z4) +fix-1-175875602+)) (setf tmp4 (mul tmp4 +fix-0-298631336+)) (setf tmp5 (mul tmp5 +fix-2-053119869+)) (setf tmp6 (mul tmp6 +fix-3-072711026+)) (setf tmp7 (mul tmp7 +fix-1-501321110+)) (setf z1 (mul z1 (- +fix-0-899976223+))) (setf z2 (mul z2 (- +fix-2-562915447+))) (setf z3 (mul z3 (- +fix-1-961570560+))) (setf z4 (mul z4 (- +fix-0-390180644+))) (incf z3 z5) (incf z4 z5) (setf (svref dptr 7) (descale (plus3 tmp4 z1 z3) +shift-1+)) (setf (svref dptr 5) (descale (plus3 tmp5 z2 z4) +shift-1+)) (setf (svref dptr 3) (descale (plus3 tmp6 z2 z3) +shift-1+)) (setf (svref dptr 1) (descale (plus3 tmp7 z1 z4) +shift-1+))) (loop for cnt fixnum from 7 downto 0 do ;second pass: on columns (setf tmp0 (plus (dbref data cnt 0) (dbref data cnt 7))) (setf tmp7 (minus (dbref data cnt 0) (dbref data cnt 7))) (setf tmp1 (plus (dbref data cnt 1) (dbref data cnt 6))) (setf tmp6 (minus (dbref data cnt 1) (dbref data cnt 6))) (setf tmp2 (plus (dbref data cnt 2) (dbref data cnt 5))) (setf tmp5 (minus (dbref data cnt 2) (dbref data cnt 5))) (setf tmp3 (plus (dbref data cnt 3) (dbref data cnt 4))) (setf tmp4 (minus (dbref data cnt 3) (dbref data cnt 4))) (setf tmp10 (plus tmp0 tmp3)) (setf tmp13 (minus tmp0 tmp3)) (setf tmp11 (plus tmp1 tmp2)) (setf tmp12 (minus tmp1 tmp2)) (setf (dbref data cnt 0) (descale (plus tmp10 tmp11) 1)) (setf (dbref data cnt 4) (descale (minus tmp10 tmp11) 1)) (setf z1 (mul (plus tmp12 tmp13) +fix-0-541196100+)) (setf (dbref data cnt 2) (descale (plus z1 (mul tmp13 +fix-0-765366865+)) +shift+1+)) (setf (dbref data cnt 6) (descale (plus z1 (mul tmp12 (- +fix-1-847759065+))) +shift+1+)) (setf z1 (plus tmp4 tmp7)) (setf z2 (plus tmp5 tmp6)) (setf z3 (plus tmp4 tmp6)) (setf z4 (plus tmp5 tmp7)) (setf z5 (mul (plus z3 z4) +fix-1-175875602+)) (setf tmp4 (mul tmp4 +fix-0-298631336+)) (setf tmp5 (mul tmp5 +fix-2-053119869+)) (setf tmp6 (mul tmp6 +fix-3-072711026+)) (setf tmp7 (mul tmp7 +fix-1-501321110+)) (setf z1 (mul z1 (- +fix-0-899976223+))) (setf z2 (mul z2 (- +fix-2-562915447+))) (setf z3 (mul z3 (- +fix-1-961570560+))) (setf z4 (mul z4 (- +fix-0-390180644+))) (incf z3 z5) (incf z4 z5) (setf (dbref data cnt 7) (descale (plus3 tmp4 z1 z3) +shift+1+)) (setf (dbref data cnt 5) (descale (plus3 tmp5 z2 z4) +shift+1+)) (setf (dbref data cnt 3) (descale (plus3 tmp6 z2 z3) +shift+1+)) (setf (dbref data cnt 1) (descale (plus3 tmp7 z1 z4) +shift+1+))) (return))) ;;; Forward DCT and quantization (defun crunch (buf pos table) (declare #.*optimize* (type fixnum pos) (type (simple-vector *) buf)) (llm-dct (svref buf pos)) (quantize (svref buf pos) table)) ;;; Q-tables scaling (defun q-scale (table q-factor) (declare #.*optimize*) (when (/= q-factor 64) (let ((factor (/ q-factor 64))) (loop for q-row across table do (loop for x fixnum from 0 to 7 do (setf (svref q-row x) (the fixnum (round (* (svref q-row x) factor))))))))) ;;; Function that maps value into SSSS (defun csize (n) (declare #.*optimize* (type fixnum n val LSB MSB)) (svref +csize+ (plus n 1023))) ;;; zigzag ordering (defun zigzag (buffer) (declare #.*optimize* (type (simple-vector 8) buffer)) (loop for row across buffer for z-row across +zigzag-index+ do (loop for x fixnum from 0 to 7 do (setf (svref *zz-result* (svref z-row x)) (the fixnum (svref row x))))) *zz-result*) ;;; Writes frame header (defun write-frame-header (maxX maxY cn q-tables sampling tqv out-stream) (declare #.*optimize* (type fixnum maxX maxY cn)) (write-huffman-tables out-stream) (write-quantization-tables q-tables out-stream) ;;writing frame header (write-marker +M_SOF0+ out-stream) (write-byte 0 out-stream) ;length (write-byte (plus 8 (mul 3 cn)) out-stream) (write-byte 8 out-stream) ;sample precision (write-byte (ash maxY -8) out-stream) ;max height (write-byte (logand maxY #xff) out-stream) (write-byte (ash maxX -8) out-stream) ;max width (write-byte (logand maxX #xff) out-stream) (write-byte cn out-stream) ;number of components (loop for entry in sampling for i fixnum from 0 by 1 do (write-byte i out-stream) (write-byte ;H and V (deposit-field (second entry) (byte 4 0)(ash (the fixnum (first entry)) 4)) out-stream) (write-byte (svref tqv i) out-stream))) ;Tq ;;; Writes byte with stuffing (adds zero after FF code) (defun write-stuffed (b s) (declare #.*optimize* (type fixnum b) (type stream s)) (write-byte b s) (if (= b #xFF) (write-byte 0 s))) ;;; A function for bit streaming ;;; NB: probably it's a good idea to encapsulate this behavior into a class, but I'm afraid that method dispatch would be too slow (defun write-bits (bi ni s) (declare #.*optimize* (special *prev-length* *prev-byte*) (type fixnum bi ni *prev-length* *prev-byte*) (type stream s)) (loop with lim fixnum = (if (> ni 8) 1 0) for i fixnum from lim downto 0 do (let ((b (ldb (byte 8 (ash i 3)) bi)) (n (cond ((and (= i 1) (= ni 16)) 8) ((and (= i 0) (/= lim 0)) 8) ((= ni 8) 8) (t (logand ni 7))))) (declare (type fixnum b n) (dynamic-extent b n)) (cond ((zerop n)) ((>= (plus n *prev-length*) 8) (let* ((result (the fixnum (ash (the fixnum *prev-byte*) (the fixnum (minus 8 *prev-length*))))) (total-length (plus n *prev-length*)) (overflow (minus total-length 8))) (declare (type fixnum overflow total-length result) (dynamic-extent overflow total-length result)) (setf *prev-byte* (ldb (byte overflow 0) b)) (write-stuffed (deposit-field (ldb (byte (minus n overflow) overflow) b) (byte (minus 8 *prev-length*) 0) result) s) (setf *prev-length* overflow))) (t (setf *prev-byte* (deposit-field b (byte n 0) (the fixnum (ash (the fixnum *prev-byte*) (the fixnum n))))) (incf *prev-length* n)))))) ;;; Encodes block using specified huffman tables, returns new pred (DC prediction value) ;;; and last code written to stream for padding (defun encode-block (block tables pred s) (declare #.*optimize* (type (simple-vector *) block)) (let* ((ehufsi-dc (first (first tables))) (ehufco-dc (second (first tables))) (ehufsi-ac (first (second tables))) (ehufco-ac (second (second tables))) (newpred (svref block 0)) (diff (minus newpred pred)) (dcpos (csize diff))) (declare (type fixnum pred newpred diff pos dcpos) (dynamic-extent diff dcpos)) ;; writing dc code first (write-bits (svref ehufco-dc dcpos) (svref ehufsi-dc dcpos) s) (cond ((minusp diff) (write-bits (1- diff) (csize diff) s)) (t (write-bits diff (csize diff) s))) ;;writing ac sequence (loop with r fixnum = 0 for k fixnum from 1 to 63 do (if (zerop (svref block k)) (if (= k 63) (progn (write-bits (svref ehufco-ac 0) (svref ehufsi-ac 0) s) ;writing EOB (return)) (incf r)) (progn (loop while (> r 15) do (write-bits (svref ehufco-ac #xf0) (svref ehufsi-ac #xf0) s) (decf r 16)) (let* ((ssss (csize (svref block k))) (rs (plus ssss (ash r 4)))) (write-bits (svref ehufco-ac rs) (svref ehufsi-ac rs) s) (when (minusp (svref block k)) (decf (svref block k) 1)) (write-bits (svref block k) ssss s)) (setf r 0)))) newpred)) ;;; Emits q-tables (defun write-quantization-tables (tables s) (let ((len (plus 2 (mul 65 (length tables))))) (write-marker +M_DQT+ s) (write-byte (ash len -8) s) ;;;MSB (write-byte (logand len #xff) s) ;;;LSB (loop for table across tables for i fixnum from 0 do (write-byte i s) (write-sequence (zigzag table) s)))) ;;; Emits huffman tables in the following order: ;;; luminance DC ;;; luminance AC ;;; chrominance DC ;;; chrominance AC (defun write-huffman-tables (s) (let ((len (+ 2 (* 17 4) (length +luminance-dc-values+) (length +luminance-ac-values+) (length +chrominance-dc-values+) (length +chrominance-ac-values+)))) (write-marker +M_DHT+ s) (write-byte (ash len -8) s) ;;;MSB (write-byte (logand len #xff) s) ;;;LSB (write-hufftable +luminance-dc-bits+ +luminance-dc-values+ 0 s) (write-hufftable +luminance-ac-bits+ +luminance-ac-values+ 16 s) (write-hufftable +chrominance-dc-bits+ +chrominance-dc-values+ 1 s) (write-hufftable +chrominance-ac-bits+ +chrominance-ac-values+ 17 s))) ;;; Writes single huffman table (defun write-hufftable (bits vals tcti s) (declare (type fixnum tcti)) (write-byte tcti s) ;Tc/Th (write-sequence bits s) (write-sequence vals s)) ;;; Drops specified marker into the stream (defun write-marker (code to) (write-sequence `(#xFF ,code) to)) ;;; Writing some markers into the stream (defun prepare-JFIF-stream (out-stream) (write-marker +M_SOI+ out-stream) (write-marker +M_APP0+ out-stream) (write-sequence '(0 16 #x4a #x46 #x49 #x46 0 1 2 0 0 1 0 1 0 0) out-stream)) ;;; Builds common decoding and encoding tables (defun build-universal-tables (bits) (let ((huffsize (make-array 256)) (huffcode (make-array 256)) (lastk 0)) (declare #.*optimize* (type fixnum lastk) (type (simple-vector *) bits huffcode huffsize)) ;; generating huffsize (loop for i fixnum from 1 to 16 with k fixnum = 0 and j fixnum = 1 do (loop until (> j (svref bits (1- i))) do (setf (svref huffsize k) i) (incf k) (incf j) finally (setf j 1)) finally (progn (setf lastk k) (setf (svref huffsize lastk) 0))) ;; generating huffcode (loop with k fixnum = 0 and code fixnum = 0 and si fixnum = (svref huffsize 0) do (loop do (setf (svref huffcode k) code) (incf code) (incf k) when (/= (svref huffsize k) si) do (return)) when (zerop (svref huffsize k)) do (return) else do (loop do (setf code (ash code 1)) (incf si) when (= (svref huffsize k) si) do (return))) (values huffcode huffsize lastk))) ;;;Builds ordered code tables for encoder (defun build-tables (bits huffval) (let ((ehufco (make-array 256)) (ehufsi (make-array 256))) (multiple-value-bind (huffcode huffsize lastk) (build-universal-tables bits) (declare (type (simple-vector *) huffsize huffcode) (type fixnum lastk)) (loop with i fixnum for k from 0 below lastk do (setf i (svref huffval k)) (setf (svref ehufco i) (svref huffcode k)) (setf (svref ehufsi i) (svref huffsize k))) (list ehufsi ehufco)))) ;;; Main encoder function (user interface) ;;; Mutilated by David Steuber to work with the fractal data (defun encode-image (filename image h w y-comp-lut u-comp-lut v-comp-lut) (declare #.*optimize* (type fixnum h w) (type (simple-array (integer 0 255) *) image)) (with-open-file (out-stream filename :direction :output :element-type 'unsigned-byte :if-exists :supersede) (let* ((ncomp 3) (q-tabs +q-tables+) (sampling '((2 2)(1 1)(1 1))) (wd (loop for entry in sampling maximize (first entry))) (ht (loop for entry in sampling maximize (second entry))) (isampling (convert-sampling sampling wd ht)) (height (ash ht 3)) (width (ash wd 3)) (YUV (make-array ncomp :initial-contents (loop for i fixnum from 0 below ncomp collecting (make-array height :initial-contents (loop for j fixnum from 0 below height collecting (make-array width)))))) (sampled-buf (make-array (mul ht wd) :initial-contents (loop for b fixnum from 0 below (mul ht wd) collecting (make-array 8 :initial-contents (loop for i fixnum from 0 to 7 collecting (make-array 8)))))) (preds (make-array 3 :initial-element 0)) (tqv #(0 1 1))) (setq *prev-byte* 0) (setq *prev-length* 0) (prepare-JFIF-stream out-stream) (write-frame-header w h 3 q-tabs sampling tqv out-stream) ;frame header ;;writing scan header (write-marker +M_SOS+ out-stream) (write-sequence '(0 12 3 0 0 1 17 2 17 0 63 0) out-stream) (let ((luminance-tabset (list (build-tables +luminance-dc-bits+ +luminance-dc-values+) (build-tables +luminance-ac-bits+ +luminance-ac-values+))) (chrominance-tabset (list (build-tables +chrominance-dc-bits+ +chrominance-dc-values+) (build-tables +chrominance-ac-bits+ +chrominance-ac-values+)))) (loop for dy fixnum from 0 below h by height do (loop for dx fixnum from 0 below w by width do (multiple-value-bind (xlim ylim) (decompose-yuv image YUV dx dy h w height width y-comp-lut u-comp-lut v-comp-lut) (declare (type fixnum xlim ylim) (dynamic-extent xlim ylim)) (loop for comp across YUV for freq in sampling for ifreq across isampling for iH fixnum = (first ifreq) for iV fixnum = (second ifreq) for cn fixnum from 0 for hufftabs = (if (zerop cn) luminance-tabset chrominance-tabset) for q-tab = (svref q-tabs (svref tqv cn)) ;choosing appropriate q-table for a component for H fixnum = (first freq) for V fixnum = (second freq) do (subsample comp sampled-buf H V (minus xlim dx) (minus ylim dy) iH iV) (loop for y fixnum from 0 below V for ypos fixnum = (if (> (plus dy (ash y 3)) ylim) (mul (rem (ash ylim -3) V) H) (mul y H)) do (loop for x fixnum from 0 below H for pos fixnum = (if (> (plus dx (ash x 3)) xlim) (plus (rem (ash xlim -3) H) ypos) (plus x ypos)) do (crunch sampled-buf pos q-tab) (setf (svref preds cn) (encode-block (zigzag (svref sampled-buf pos)) hufftabs (svref preds cn) out-stream))))))))) (unless (zerop *prev-length*) (write-stuffed (deposit-field #xff ;byte padding & flushing (byte (minus 8 *prev-length*) 0) (the fixnum (ash (the fixnum *prev-byte*) (the fixnum (minus 8 *prev-length*))))) out-stream)) (write-marker +M_EOI+ out-stream)))) ;;;; ************************************************************************ ;;;; ************************************************************************ ;;;; ************************************************************************ ;;;; ************************************************************************ ;;;; ************************************************************************ ;;;; Mandelbrot Code (in-package :cl-user) (eval-when (:compile-toplevel :load-toplevel :execute) (defconstant +max-iterations+ 128)) (defclass mset-map () ((left :type single-float :initarg :left :reader mset-map-left) (right :type single-float :initarg :right :reader mset-map-right) (top :type single-float :initarg :top :reader mset-map-top) (bottom :type single-float :initarg :bottom :reader mset-map-bottom) (width :type (integer 1 4096) :initarg :width :reader mset-map-width) (height :type (integer 1 4096) :initarg :height :reader mset-map-height) (x-axis :type simple-vector) (y-axis :type simple-vector) (data :type (array (integer 0 #.+max-iterations+) (* *))) (data-alias :type (simple-array (integer 0 #.+max-iterations+) *)))) (defmethod print-object ((m mset-map) s) (with-slots (left right top bottom width height) m (print-unreadable-object (m s :type t :identity t) (format s "[~A ~A ~A ~A] ~Ax~A" left right top bottom width height)))) (defmethod shared-initialize :after ((o mset-map) slot-names &rest init-args) (declare (ignore slot-names init-args)) (flet ((make-axis-array (len min max) (let ((my-array (make-array len))) (loop for i below len with pixelsize = (/ (- max min) len) with start = min do (setf (svref my-array i) (+ start (* i pixelsize)))) my-array))) (with-slots (left right top bottom x-axis y-axis data data-alias height width) o (setf x-axis (make-axis-array width left right) y-axis (nreverse (make-axis-array height bottom top)) data-alias (make-array (* height width) :element-type '(integer 0 #.+max-iterations+) :initial-element 0) data (make-array (list height width) :element-type '(integer 0 #.+max-iterations+) :displaced-to data-alias))))) (defun mset-map-clear-data (mset-map) (let ((data-alias (slot-value mset-map 'data-alias))) (declare (type (simple-array (integer 0 #.+max-iterations+) *) data-alias) (optimize speed)) (loop for i of-type (unsigned-byte 24) below (length data-alias) do (setf (aref data-alias i) 0)))) (defun mset-map-reset-coords (mset-map left right bottom top) (declare (type single-float left right bottom top) (optimize speed)) (flet ((reset-axis-array (array len min max) (declare (type (simple-vector *) array) (type (integer 1 1024) len) (type single-float min max)) (loop with pixelsize of-type single-float = (/ (- max min) len) with start of-type single-float = min for i of-type (integer 1 1024) below len do (setf (svref array i) (+ start (* i pixelsize)))) array) (single-float (n) (coerce n 'single-float))) (mset-map-clear-data mset-map) (with-slots (x-axis y-axis width height (l left) (r right) (to top) (b bottom)) mset-map (setf l (single-float left) r (single-float right) to (single-float top) b (single-float bottom)) (setf x-axis (reset-axis-array x-axis width l r) y-axis (nreverse (reset-axis-array y-axis height b to))))) mset-map) (defun escape-fn (c-real c-imag) "Return iterations for c-real + c-imag i as floating-point numbers" (declare (type single-float c-real c-imag) (optimize (speed 3) (safety 0) (debug 0))) (loop for iterations of-type (integer 0 #.+max-iterations+) below +max-iterations+ for period of-type (integer 0 #.(ash +max-iterations+ 1)) from 1 for z-imag of-type single-float = 0.0 then (+ (* 2.0 z-real z-imag) c-imag) for z-real of-type single-float = 0.0 then (+ (- z-real-squared z-imag-squared) c-real) for z-real-squared of-type single-float = 0.0 then (* z-real z-real) for z-imag-squared of-type single-float = 0.0 then (* z-imag z-imag) with a-imag of-type single-float = 0.0 with a-real of-type single-float = 0.0 with check of-type (integer 0 #.+max-iterations+) = 3 unless (< (+ z-real-squared z-imag-squared) 4.0) return iterations when (and (= a-imag z-imag) (= a-real z-real) (> iterations 1)) return +max-iterations+ when (> period check) do (setf period 1 check (ash check 1) a-imag z-imag a-real z-real) finally (return iterations))) (defgeneric compute-row (mset-map row left right)) (defmethod compute-row ((m mset-map) row left right) (declare (type (unsigned-byte 16) row left right) (optimize speed)) (let ((x-axis (slot-value m 'x-axis)) (y-axis (slot-value m 'y-axis)) (data (slot-value m 'data))) (declare (type (simple-array (integer 0 #.+max-iterations+) (* *)) data) (type (simple-vector *) x-axis y-axis)) (loop with y of-type single-float = (the single-float (svref y-axis row)) for col from left to right when (= 0 (the (integer 0 #.+max-iterations+) (aref data row col))) do (setf (aref data row col) (the (integer 0 #.+max-iterations+) (escape-fn (svref x-axis col) y)))))) (defgeneric compute-column (mset-map col top bottom)) (defmethod compute-column ((m mset-map) col top bottom) (declare (type (unsigned-byte 16) col top bottom) (optimize speed)) (let ((x-axis (slot-value m 'x-axis)) (y-axis (slot-value m 'y-axis)) (data (slot-value m 'data))) (declare (type (simple-array (integer 0 #.+max-iterations+) (* *)) data) (type (simple-vector *) x-axis y-axis)) (loop with x of-type single-float = (the single-float (svref x-axis col)) for row from top to bottom when (= 0 (the (integer 0 #.+max-iterations+) (aref data row col))) do (setf (aref data row col) (the (integer 0 #.+max-iterations+) (escape-fn x (svref y-axis row))))))) (defgeneric check-edges (mset-map left top right bottom) (:documentation "If all the values are the same, fill the region and return t")) (defmethod check-edges ((m mset-map) left top right bottom) (declare (type (unsigned-byte 16) left top right bottom) (optimize (speed 3) (safety 0) (debug 0))) (let ((data (slot-value m 'data))) (declare (type (simple-array (integer 0 #.+max-iterations+) (* *)) data)) (flet ((fill-region (v) (declare (type (integer 0 #.+max-iterations+) v)) (loop for row of-type (unsigned-byte 16) from (1+ top) to (1- bottom) do (loop for col of-type (unsigned-byte 16) from (1+ left) to (1- right) do (setf (aref data row col) v))))) (let ((v (the (integer 0 #.+max-iterations+) (aref data top left))) (r t)) (declare (type (integer 0 #.+max-iterations+) v)) (loop for col from left to right unless (= v (aref data top col) (aref data bottom col)) return (setf r nil)) (when r (loop for row from top to bottom unless (= v (aref data row left) (aref data row right)) return (setf r nil))) (when r (fill-region v)) r)))) (defgeneric quarter-the-map (mset-map left top right bottom) (:documentation "Divide the map into four quarters")) (defmethod quarter-the-map ((m mset-map) left top right bottom) (declare (type (unsigned-byte 16) left top right bottom) (optimize speed)) (cond ((< (- right left) 2) nil) ; done ((< (- bottom top) 2) nil) ; also done ((check-edges m left top right bottom) nil) ; done here also (t (let ((vmid (+ top (the (unsigned-byte 16) (round (- bottom top) 2)))) (hmid (+ left (the (unsigned-byte 16) (round (- right left) 2))))) (compute-row m vmid (1+ left) (1- right)) (compute-column m hmid (1+ top) (1- bottom)) (quarter-the-map m left top hmid vmid) (quarter-the-map m hmid top right vmid) (quarter-the-map m hmid vmid right bottom) (quarter-the-map m left vmid hmid bottom))))) (defgeneric compute-mandelbrot-set (mset-map) (:documentation "Calculate bitmap data for the Mandelbrot Set")) (defmethod compute-mandelbrot-set ((m mset-map)) (with-slots (width height) m (compute-row m 0 0 (1- width)) (compute-row m (1- height) 0 (1- width)) (compute-column m 0 1 (- height 2)) (compute-column m (1- width) 1 (- height 2)) (quarter-the-map m 0 0 (1- width) (1- height))) m) (defgeneric mset-map-row-col (mset-map row col)) (defmethod mset-map-row-col ((m mset-map) row col) (aref (slot-value m 'data) row col)) ;; Color Map (defclass color-rgb () ((red :type (unsigned-byte 8) :initarg :red :initform 0 :accessor color-rgb-red) (green :type (unisgned-byte 8) :initarg :green :initform 0 :accessor color-rgb-green) (blue :type (unsigned-byte 8) :initarg :blue :initform 0 :accessor color-rgb-blue))) (defclass color-yuv (color-rgb) ()) (defgeneric color-yuv-y (color-yuv)) (defgeneric color-yuv-u (color-yuv)) (defgeneric color-yuv-v (color-yuv)) (defmethod color-yuv-y ((c color-yuv)) (with-slots (red green blue) c (jpeg::y-comp red green blue))) (defmethod color-yuv-u ((c color-yuv)) (with-slots (red green blue) c (jpeg::u-comp red green blue))) (defmethod color-yuv-v ((c color-yuv)) (with-slots (red green blue) c (jpeg::v-comp red green blue))) (defun make-color-map () (let ((cm (make-array (1+ +max-iterations+)))) (flet ((sincos-to-byte (n) (round (+ 255 (* 255 n)) 2)) (itts-to-rad (i) (* (/ (if (zerop i) 0 (log i 1.1)) +max-iterations+) 4.5 pi))) (loop for i below +max-iterations+ as c = (make-instance 'color-yuv :red (sincos-to-byte (sin (- (* 0.50 pi) (itts-to-rad i)))) :green (sincos-to-byte (cos (+ (* 0.50 pi) (itts-to-rad i)))) :blue (sincos-to-byte (sin (+ pi (itts-to-rad i))))) do (setf (aref cm i) c) finally (setf (aref cm i) (make-instance 'color-yuv :red 0 :green 0 :blue 0))) cm))) (defparameter *color-map* (make-color-map)) (defparameter *color-map-y* (make-array (1+ +max-iterations+) :element-type 'signed-byte :initial-contents (loop for y across *color-map* collect (color-yuv-y y)))) (defparameter *color-map-u* (make-array (1+ +max-iterations+) :element-type 'signed-byte :initial-contents (loop for u across *color-map* collect (color-yuv-u u)))) (defparameter *color-map-v* (make-array (1+ +max-iterations+) :element-type 'signed-byte :initial-contents (loop for v across *color-map* collect (color-yuv-v v)))) ;; File I/O (defgeneric write-mset-map-to-jpeg-file (mset-map filespec)) (defmethod write-mset-map-to-jpeg-file ((m mset-map) filespec) (with-slots (data-alias height width) m (jpeg:encode-image filespec data-alias height width *color-map-y* *color-map-u* *color-map-v*))) (defun generate-fractals-from-input-file (filespec &optional (outdir "outputpic") (width 800) (height 600)) (let ((*read-eval* nil) (mset-map (make-instance 'mset-map :width width :height height :left -1.0 :right 1.0 :top 1.0 :bottom -1.0)) (outdir (make-pathname :directory `(:relative ,outdir)))) (ensure-directories-exist outdir) (with-open-file (s filespec :direction :input) (loop for image from 1 as outfile = (merge-pathnames (make-pathname :name (format nil "image~2,'0d" image) :type "jpeg") outdir) as coords = (loop for i below 4 as coord = (read s nil nil) while coord collect coord) while coords do (destructuring-bind (left right bottom top) coords (mset-map-reset-coords mset-map left right bottom top) (format t "File: ~S~%" outfile) (format t "left: ~A, right: ~A, bottom: ~A, top: ~A~%" (mset-map-left mset-map) (mset-map-right mset-map) (mset-map-bottom mset-map) (mset-map-top mset-map)) (compute-mandelbrot-set mset-map) (write-mset-map-to-jpeg-file mset-map outfile))))))