;;;; -*- Mode: LISP; Syntax: ANSI-Common-Lisp; Base: 10; Package: photo -*-
;;;; *************************************************************************
;;;; FILE IDENTIFICATION
;;;;
;;;; Name:          dof.lisp
;;;; Purpose:       Depth of field functions for cl-photo
;;;; Programmer:    Kevin M. Rosenberg
;;;; Date Started:  April 2005
;;;;
;;;; $Id$
;;;;
;;;; This file, part of cl-photo, is Copyright (c) 2005 by Kevin M. Rosenberg
;;;;
;;;; cl-photo users are granted the rights to distribute and use this software
;;;; as governed by the terms of the GNU General Public License v2
;;;; (http://www.gnu.org/licenses/gpl.html)
;;;;
;;;; *************************************************************************

(in-package #:photo)

(defun sort-size (size)
  "Returns a cons pair with the smaller size first."
  (if (>= (car size) (cdr size))
      (cons (cdr size) (car size))
      (cons (car size) (cdr size))))

(defun print-magnification (imager-size print-size)
  "Returns the magnification required between an imager and print sizes
while taking crop into consideration."
  (setf imager-size (sort-size imager-size))
  (setf print-size (sort-size print-size))
  (float (max (/ (car print-size) (car imager-size))
              (/ (cdr print-size) (cdr print-size)))))

(defun coc (imager-size &key (lpm 5) (minimum-distance 250)
                   (viewing-distance 250)
                   (print-size (output-dimensions :8x10in)))
  "Returns circle of confusion in mm and print magnification for a format.
Default resolving power is 5 lpm at 25cm."
  (let* ((magnification (print-magnification imager-size print-size))
         (resolution-factor (/ (* magnification lpm minimum-distance) viewing-distance))
         (coc (/ 1.0d0 resolution-factor)))
    (values coc magnification)))

(defun coc-format (format &key (lpm 5) (minimum-distance 250)
                          (viewing-distance 250)
                          (print-size (output-dimensions :8x10in)))
  "Returns circle of confusion in mm and print magnification for a format.
Default resolving power is 5 lpm at 25cm."

  (let* ((format-size (imager-dimensions format))
         (format-diagonal (diagonal (car format-size) (cdr format-size)))
         (print-diagonal (diagonal (car print-size) (cdr print-size)))
         (resolution-factor (/ (* lpm print-diagonal minimum-distance)
                               (* format-diagonal viewing-distance)))
         (coc (/ 1.0d0 resolution-factor))
         (print-magnification (/ print-diagonal format-diagonal)))
    (values coc print-magnification)))

(defun coc-pixels (imager pixels)
  "Returns lpm and circle of confusion based on pixel size."
  (when (and (consp imager) (consp pixels))
    (let ((coc-w (float (* 2 (/ (car imager) (car pixels)))))
          (coc-h (float (* 2 (/ (cdr imager) (cdr pixels))))))
    (values coc-w coc-h (/ 1. coc-w) (/ 1. coc-h)))))

(defun coc-pixels-format (format)
  "Returns circle of confusion based on pixel size."
  (coc-pixels (imager-dimensions format) (pixel-dimensions format)))

(defun coc-airy (f-stop &optional (wavelength 0.000512))
  "Return the circle of confusion based on the airy disk."
  (float (/ 1 (rayleigh-limit f-stop wavelength))))

(defun rayleigh-limit (f-stop &optional (wavelength 0.0005))
  "Returns the rayleigh limit in line pairs per mm (MTF 9%) as well as the MTF50"
  (let ((rayleigh (float (/ 1 1.22 f-stop wavelength))))
    (values rayleigh (* 0.46 rayleigh))))

(defun maximum-sharpness-aperture (format &optional (wavelength 0.0005))
  (multiple-value-bind (coc-w coc-h lpm-w lpm-h) (coc-pixels-format format)
    (declare (ignore coc-w coc-h))
    (/ 1. (* 1.22 wavelength (/ (min lpm-w lpm-h) 0.46)))))

(defun dof-mm (focal-length f-stop distance coc &key (pupil-factor 1))
  "Returns depth of field based on focal-length, f-stop, distance, and coc.
Six values are returned:
near point, far point, total dof, magnification, blur size at infinity (mm).
Circle of confusion can either be a number or keyword designating format.
Reference: http://www.vanwalree.com/optics/dofderivation.html"
  (let* ((aperture (/ focal-length f-stop))
         (hyperfocal (hyperfocal focal-length f-stop coc))
         (numerator-1 (* (- pupil-factor 1) (- distance focal-length)
                         coc focal-length))
         (numerator-2 (* pupil-factor aperture focal-length distance))
         (denominator-1 (* pupil-factor coc (- distance focal-length)))
         (denominator-2 (* pupil-factor aperture focal-length))
         (near (/ (+ numerator-1 numerator-2)
                  (+ denominator-1 denominator-2)))
         (far (when (/= denominator-1 denominator-2)
                (/ (- numerator-1 numerator-2)
                   (- denominator-1 denominator-2))))
         (mag (float (/ focal-length (- distance focal-length))))
         (infinity-blur-diameter (/ (* mag focal-length) f-stop))
         (depth (when far (- far near))))
    (when (or (>= distance hyperfocal)
               (and (null far) (>= distance (* hyperfocal 0.99))))
      (setq near (/ hyperfocal 2)
            far most-positive-short-float
            depth most-positive-short-float))
    (values near far depth mag infinity-blur-diameter)))

;; Simplified calculation for symmetric lens
(defun dof-symmetric-mm (focal-length f-stop distance coc)
  "Returns depth of field based on focal-length, f-stop, distance, and coc.
Six values are returned:
near point, far point, total dof, near point, far point, magnification,
blur size at infinity (mm).
Circle of confusion can either be a number or keyword designating format."
  (let* ((aperture (/ focal-length f-stop))
         (hyperfocal (hyperfocal focal-length f-stop coc))
         (numerator (* distance coc (- distance focal-length)))
         (factor-1 (* focal-length aperture))
         (factor-2 (* coc (- distance focal-length)))
         (near (- distance (/ numerator (+ factor-1 factor-2))))
         (far (when (/= factor-1 factor-2)
                (+ distance (/ numerator (- factor-1 factor-2)))))
         (mag (magnification :focal-length focal-length :object-distance distance :units :mm))
         (infinity-blur-diameter (/ (* mag focal-length) f-stop))
         (depth (when far (- far near))))
    (when (or (>= distance hyperfocal)
               (and (null far) (>= distance (* hyperfocal 0.99))))
      (setq near (/ hyperfocal 2)
            far most-positive-short-float
            depth most-positive-short-float))
    (values near far depth mag infinity-blur-diameter)))

(defun dof (focal-length f-stop distance coc &key (units :mm) (pupil-factor 1))
  "Returns the Depth of Field.
Input: FOCAL-LENGTH, F-STOP, DISTANCE, CIRCLE-OF-CONFUSION.
Output: NEAR-POINT, FAR-POINT, TOTAL-DOF, MAGNIFICATION, BLUR-SIZE-OF-INFINITY-POINT-IN-MM."
  (multiple-value-bind (near-point far-point total-dof mag blur)
      (dof-mm focal-length f-stop (length->mm distance units) coc
              :pupil-factor pupil-factor)
    (values (mm->length near-point units)
            (mm->length far-point units)
            (mm->length total-dof units)
            mag blur)))

(defun hyperfocal (focal-length f-stop coc &key (units :mm))
  (mm->length (+ focal-length (/ (* focal-length focal-length) f-stop coc)) units))

(defun effective-aperture (focal-length distance aperture)
  (* aperture (bellows-factor focal-length distance)))

(defun mtf-scanner (freq dscan-freq &optional (order 3))
  (abs (expt (kmrcl:sinc (* pi (/ freq dscan-freq))) order)))