home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / ascendMar8.tar / UMass / LineFinder / topdown-line-finder.lisp < prev   
Lisp/Scheme  |  1996-03-05  |  16KB  |  432 lines
  1. ;;; topdown-line-finder.lisp
  2. ;;; find lines with a particular orientation, in an image subarea
  3. ;;;
  4. ;;; relies on routines in RCDE for computing canny edges
  5. ;;; relies on routines in ISR2 for computing connected components
  6. ;;;
  7. ;;; Author: Bob Collins  Date: 7/5/95
  8. ;;; Copyright: University of Massachusetts, 1995, all rights reserved
  9. ;;;
  10. ;;; Modifications:
  11. ;;;  Bob Collins Tue Nov  7 14:43:31 EST 1995
  12. ;;;    Changed transform-2d-bounding-box to work correctly for transforms 
  13. ;;;    that rotate the image (actually, I have changed it so it will work 
  14. ;;;    form arbitary transforms).  It now explicitly transforms all four 
  15. ;;;    corners of the old bounding box, then computes a new bounding box 
  16. ;;;    that circumscribes this quadrilateral.  Added function
  17. ;;;    transform-2d-angle for the same reason.
  18.  
  19. (in-package 'cme)
  20.  
  21. (defstruct edgeimg
  22.   (image nil)      ;subimage for which edges were computed
  23.   (bit nil)        ;binary image containing on/off edges
  24.   (gradx nil)      ;x component of gradient
  25.   (grady nil)      ;y component of gradient
  26.   (mag nil)        ;magnitude of gradient
  27.   (filtered nil))  ;temp space for filtered canny edgels
  28.  
  29.  
  30. ;;;=============================================================================
  31. ;;; CANNY EDGE DETECTION
  32. ;;;the following is a rewrite of canny-edges from $CMEHOME/ic/canny.lisp 
  33. ;;;in order to return the images I want.
  34.  
  35. (in-package 'canny)
  36.  
  37. (defvar *mask-size* 4 "size of canny mask")
  38. (defvar *low-threshold* 1.0 "canny low threshold")
  39. (defvar *high-threshold* 4.0 "canny high threshold")
  40.  
  41. (defun cme::extract-canny-edges-and-grads (image &key
  42.                        (mask-size *mask-size*)
  43.                        (lo-threshold *low-threshold*)
  44.                        (hi-threshold *high-threshold*))
  45.   "  Returns in a structure the canny edges, the x components of 
  46.   the gradient vector, the y components of the gradient vector, 
  47.   and the gradient magnitudes."
  48.   (multiple-value-bind (x-image y-image mag-image smooth)
  49.       (grads-phase image mask-size)
  50.       (ic::unmake-image smooth)
  51.     (let ((edges (nms-phase x-image y-image mag-image
  52.                 lo-threshold hi-threshold)))
  53.       (cme::make-edgeimg
  54.        :image image
  55.        :bit edges
  56.        :gradx x-image
  57.        :grady y-image
  58.        :mag mag-image
  59.        :filtered nil))))
  60.  
  61.  
  62. (in-package 'cme)
  63.  
  64. ;;;=============================================================================
  65. ;;; FREE UP IMAGE MEMORY
  66.  
  67. (defun free-edgeimg (edges)
  68.   (when edges
  69.     (when (edgeimg-image edges) (ic::unmake-image (edgeimg-image edges)))
  70.     (when (edgeimg-bit edges) (ic::unmake-image (edgeimg-bit edges)))
  71.     (when (edgeimg-gradx edges) (ic::unmake-image (edgeimg-gradx edges)))
  72.     (when (edgeimg-grady edges) (ic::unmake-image (edgeimg-grady edges)))
  73.     (when (edgeimg-mag edges) (ic::unmake-image (edgeimg-mag edges)))
  74.     (when (edgeimg-filtered edges) (ic::unmake-image (edgeimg-filtered edges))))
  75.   edges)
  76.  
  77. ;;;=============================================================================
  78. ;;; THRESHOLDING EDGES BASED ON ORIENTATION
  79.  
  80. (defun filter-edges-by-orientation (edges mean-grad-angle
  81.             &key (delta-angle 20) (directional? nil)
  82.             (startx 0)(starty 0)
  83.             (endx (1- (image-x-dim (edgeimg-bit edges))))
  84.             (endy (1- (image-y-dim (edgeimg-bit edges)))))
  85.   "Returns a new binary edge image that only contains edges passing an orientation 
  86.   test, namely they gradient at that pixel must fall within delta-angle degrees of
  87.   mean-grad-angle.  Normally we filter on unsigned orientations, so that all
  88.   angles are considered to be (mod angle 180), however if directional-p? is
  89.   nonnil, all angles will be considered to be (mod angle 360).  Input is a 
  90.   structure of type edgeimg, output is a binary image."
  91.   (let ((edgebits (edgeimg-bit edges))
  92.     (gradx (edgeimg-gradx edges))
  93.     (grady (edgeimg-grady edges))
  94.     (gradmag (edgeimg-mag edges))
  95.     (filtered-edges (edgeimg-filtered edges)))
  96.     (unless filtered-edges
  97.        (setf filtered-edges (ic::similar-image edgebits))
  98.        (setf (edgeimg-filtered edges) filtered-edges))
  99.     (setf startx (max 0 startx) starty (max 0 starty)
  100.       endx (min (1- (image-x-dim edgebits)) endx)
  101.       endy (min (1- (image-y-dim edgebits)) endy))
  102.     (flet ((deg2rad (deg) (/ (* deg pi) 180.0)))
  103.       (setf mean-grad-angle (deg2rad mean-grad-angle))
  104.       (setf delta-angle (deg2rad delta-angle))
  105.       (let ((unit-mean-x (cos mean-grad-angle))
  106.         (unit-mean-y (sin mean-grad-angle))
  107.         (cosdelta (cos (abs delta-angle))))
  108.     (do ((y starty (1+ y)))
  109.         ((> y endy) nil)
  110.         (do ((x startx (1+ x)))
  111.         ((> x endx) nil)
  112.           (setf (iref filtered-edges x y)
  113.          (if (= (iref edgebits x y) 1)
  114.              (let* ((mag (iref gradmag x y))
  115.                 (unit-x (/ (iref gradx x y) mag))
  116.                 (unit-y (/ (iref grady x y) mag))
  117.                 (cosangle (+ (* unit-x unit-mean-x)
  118.                      (* unit-y unit-mean-y))))
  119.                (if (< (if directional? cosangle (abs cosangle))
  120.                   cosdelta)
  121.                0
  122.                1))
  123.            0))))))
  124.     filtered-edges))
  125.  
  126. ;;;=========================================================================
  127. ;;; CONNECTED COMPONENTS
  128.  
  129. ;;; the following functions use code in $CME/code/UMass/ISR/
  130. ;;; for generating line adjacency graphs.
  131.  
  132. (defun construct-edgel-lag (edgeimage startx starty endx endy 8-connect?)
  133.   "Construct a line adjacency graph of the edgels."
  134.   (let ((lag (isr2::make-lag 0 (- endy starty -1))))
  135.     (do ((y starty (1+ y))
  136.      (in-segment nil nil)
  137.      (segment-start))
  138.     ((> y endy) nil)
  139.      (do ((x startx (1+ x)))
  140.      ((> x endx) nil)
  141.        (cond
  142.      ((and in-segment (zerop (iref edgeimage x y)))
  143.         (if 8-connect?
  144.         (isr2::enter-8-connected-segment-into-lag
  145.            lag (- y starty) segment-start (1- x))
  146.             (isr2::enter-4-connected-segment-into-lag
  147.                lag (- y starty) segment-start (1- x)))
  148.         (setf in-segment nil))
  149.      ((and (not in-segment) (= (iref edgeimage x y) 1))
  150.         (setf segment-start x)
  151.         (setf in-segment t))
  152.      (t nil)))
  153.       (when in-segment
  154.      (if 8-connect?
  155.          (isr2::enter-8-connected-segment-into-lag
  156.             lag (- y starty) segment-start endx)
  157.          (isr2::enter-4-connected-segment-into-lag
  158.             lag (- y starty) segment-start endx))))
  159.     lag))
  160.  
  161.  
  162. (defun extract-edgel-chains (edgeimage &key (8-connect? t) (size-thresh 3)
  163.                   (startx 0)(starty 0)
  164.                   (endx (1- (image-x-dim edgeimage)))
  165.                   (endy (1- (image-y-dim edgeimage))))
  166.   "Extract edgel chains from edge image via connected components."
  167.   (setf startx (max 0 startx) starty (max 0 starty)
  168.     endx (min (1- (image-x-dim edgeimage)) endx)
  169.     endy (min (1- (image-y-dim edgeimage)) endy))
  170.   (unless (or (> startx endx) (> starty endy))
  171.     (let ((lag (construct-edgel-lag edgeimage 
  172.            startx starty endx endy 8-connect?))
  173.     (list-of-chains nil)
  174.     (pixlist nil))
  175.     (isr2::clear-lag-segment-marks lag)
  176.     (isr2::for-each-unmarked-lag-segment (lag ycoord segment)
  177.        (setf pixlist nil)
  178.        (isr2::mark-each-connected-segment
  179.       segment ycoord
  180.       :function
  181.          #'(lambda (segment y) 
  182.          (do ((x (isr2::segment-entry-startx segment) (1+ x)))
  183.              ((> x (isr2::segment-entry-endx segment)))
  184.            (push (list x (+ y starty)) pixlist))))
  185.        (when (> (length pixlist) size-thresh)
  186.          (push pixlist list-of-chains)))
  187.     list-of-chains)))
  188.  
  189. ;;;=========================================================================
  190. ;;; LINE FITTING TO A SET OF EDGELS -- SIMPLE LEAST-SQUARES
  191.  
  192. (defun max-unit-eigenvector (a b c)
  193.   "  Returns x and y components of the unit eigenvector associated
  194.   with the largest eigenvalue of symmetric matrix ((a b)(b c))"
  195.   (if (zerop b)
  196.       (if (<= a c) 
  197.       (values 0.0 1.0)
  198.       (values 1.0 0.0))
  199.     ;;the following calculation is based on an article by T. A. Newton in 
  200.     ;;the Jan 1990 issue of the American Mathematical Monthly.
  201.     (let* ((me (/ (- c a) b))
  202.        (m (if (minusp b)
  203.           (/ (+ me (sqrt (+ (* me me) 4.0))) 2.0)
  204.           (/ (- me (sqrt (+ (* me me) 4.0))) 2.0)))
  205.        (mag (sqrt (+ (* m m) 1.0))))
  206.       (values (/ m mag) (- (/ 1.0 mag))))))
  207.  
  208.  
  209. #|
  210. (defun eigen-2by2 (a b c)
  211.   "  Eigendecomposition of a 2 by 2 symmetric matrix ((a b)(b c)).  
  212.   Multiple values returned are max eigenvalue, min eigenvalue, 
  213.   x component of max unit eigenvector, y component of max unit eigenvector."
  214.   (let* ((trace (+ a c))
  215.          (det (- (* a c) (* b b)))
  216.          (temp (sqrt (- (* trace trace) (* 4 det))))
  217.          (eigmax (/ (+ trace temp) 2))
  218.          (eigmin (/ (- trace temp) 2))
  219.          (dely (- eigmax a))
  220.          (delx b)
  221.      (norm (sqrt (+ (* delx delx) (* dely dely)))))
  222.     (if (zerop delx)
  223.     (values eigmax eigmin 1 0)
  224.         (values eigmax eigmin (/ delx norm) (/ dely norm)))))
  225. |#
  226.  
  227. (defun fit-line-to-pixlist (pixlist &key (pixextension 1.0))
  228.   "Least squares fitting of a line segment to a list of pixels coordinates."
  229.   (let ((xsum 0)(ysum 0)
  230.     (numsamples (length pixlist))) 
  231.     (dolist (xypair pixlist)
  232.        (incf xsum (first xypair))
  233.        (incf ysum (second xypair)))
  234.     (let ((xmean (/ xsum (float numsamples)))
  235.       (ymean (/ ysum (float numsamples)))
  236.       (xx 0.0)(xy 0.0)(yy 0.0))
  237.       (dolist (xypair pixlist)
  238.      (let ((x (- (first xypair) xmean))
  239.            (y (- (second xypair) ymean)))
  240.        (incf xx (* x x))
  241.        (incf xy (* x y))
  242.        (incf yy (* y y))))
  243.       (multiple-value-bind (ux uy) (max-unit-eigenvector xx xy yy)
  244.     (let ((minproj 999999.9)(maxproj -999999.9))
  245.       (mapcar #'(lambda (xy)
  246.               (let ((proj (+ (* ux (- (car xy) xmean))
  247.                      (* uy (- (cadr xy) ymean)))))
  248.             (when (< proj minproj) (setf minproj proj))
  249.             (when (> proj maxproj) (setf maxproj proj))))
  250.           pixlist)
  251.       (incf minproj (* (signum minproj) pixextension))
  252.        (incf maxproj (* (signum maxproj) pixextension))
  253.       (list (+ xmean (* minproj ux) 0.5)
  254.         (+ ymean (* minproj uy) 0.5)
  255.         (+ xmean (* maxproj ux) 0.5)
  256.         (+ ymean (* maxproj uy) 0.5)))))))
  257.  
  258. (defun transform-linelist-coordinates (linelist transform)
  259.   (mapcar #'(lambda (endpts)
  260.           (bind-vector-elements 
  261.              (u1 v1)(isr2rcde-inline-transform 
  262.                transform (list (car endpts) (cadr endpts)))
  263.          (bind-vector-elements
  264.             (u2 v2)(isr2rcde-inline-transform 
  265.                   transform (list (third endpts) (fourth endpts)))
  266.             (list u1 v1 u2 v2))))
  267.         linelist))
  268.   
  269. ;;;=========================================================================
  270. ;;; FIND ORIENTED LINE SEGMENT FEATURES
  271.  
  272. (defun transform-2d-bounding-box (transform x1 y1 x2 y2)
  273.   (bind-vector-elements 
  274.      (u1 v1)(isr2rcde-inline-transform transform (list x1 y1))
  275.      (bind-vector-elements 
  276.         (u2 v2)(isr2rcde-inline-transform transform (list x2 y2))
  277.     (bind-vector-elements 
  278.            (u3 v3)(isr2rcde-inline-transform transform (list x2 y1))
  279.        (bind-vector-elements 
  280.               (u4 v4)(isr2rcde-inline-transform transform (list x1 y2))
  281.           (values (min u1 u2 u3 u4)
  282.               (min v1 v2 v3 v4)
  283.               (max u1 u2 u3 u4)
  284.               (max v1 v2 v3 v4)))))))
  285.  
  286. (defun transform-2d-angle (transform x y angle-degrees)
  287.   (let ((rad (radians angle-degrees)))
  288.     (let ((x1 x)
  289.       (y1 y)
  290.       (x2 (+ x (cos rad)))
  291.       (y2 (+ y (sin rad))))
  292.       (bind-vector-elements 
  293.         (u1 v1)(isr2rcde-inline-transform transform (list x1 y1))
  294.     (bind-vector-elements 
  295.           (u2 v2)(isr2rcde-inline-transform transform (list x2 y2))
  296.       (degrees (atan (- v2 v1) (- u2 u1))))))))
  297.  
  298. (defun transform-2dworld-bounding-box (image startx starty endx endy)
  299.   (let ((transform (inverse-transform (image-to-2d-transform image))))
  300.     (multiple-value-bind (u1 v1 u2 v2) 
  301.        (transform-2d-bounding-box transform startx starty endx endy)
  302.        (values (floor u1) (floor v1) (ceiling u2) (ceiling v2)))))
  303.  
  304. (defun transform-2dworld-angle (image x y angle-degrees)
  305.   (let ((transform (inverse-transform (image-to-2d-transform image))))
  306.     (transform-2d-angle transform x y angle-degrees)))
  307.  
  308.  
  309. (defun compute-2dworld-canny-edges (image startx starty endx endy
  310.                      &key (mask-size canny::*mask-size*)
  311.                       (lo-threshold canny::*low-threshold*)
  312.                       (hi-threshold canny::*high-threshold*))
  313.   "  Compute canny edges for the given image within bounding box  
  314.    startx starty endx endy, specified in 2d-world coordinates."
  315.   (let ((x1 (min startx endx))
  316.     (x2 (max startx endx))
  317.     (y1 (min starty endy))
  318.     (y2 (max starty endy)))
  319.     (multiple-value-bind (u1 v1 u2 v2) 
  320.           (transform-2dworld-bounding-box image x1 y1 x2 y2) 
  321.      (unless (or (< u2 0) (< v2 0)
  322.          (>= u1 (image-x-dim image))
  323.          (>= v1 (image-y-dim image)))
  324.       (let* ((halfsize (1+ (ash mask-size -1)))
  325.         (subimage (ic::image-window image 
  326.               (- u1 halfsize) (- v1 halfsize)
  327.               (+ (- u2 u1 -1) (* 2 halfsize) 1)
  328.               (+ (- v2 v1 -1) (* 2 halfsize) 1))))
  329.     (extract-canny-edges-and-grads 
  330.        subimage 
  331.        :mask-size mask-size
  332.        :lo-threshold lo-threshold
  333.        :hi-threshold hi-threshold))))))
  334.  
  335. (defun compute-imageuv-canny-edges (image startu startv endu endv
  336.                      &key (mask-size canny::*mask-size*)
  337.                       (lo-threshold canny::*low-threshold*)
  338.                       (hi-threshold canny::*high-threshold*))
  339.   "  Compute canny edges for the given image within bounding box  
  340.    startx starty endx endy, specified in 2d-world coordinates."
  341.   (let ((u1 (min startu endu))
  342.     (u2 (max startu endu))
  343.     (v1 (min startv endv))
  344.     (v2 (max startv endv)))
  345.      (unless (or (< u2 0) (< v2 0)
  346.          (>= u1 (image-x-dim image))
  347.          (>= v1 (image-y-dim image)))
  348.       (let* ((halfsize (1+ (ash mask-size -1)))
  349.         (subimage (ic::image-window image 
  350.               (- u1 halfsize) (- v1 halfsize)
  351.               (+ (- u2 u1 -1) (* 2 halfsize) 1)
  352.               (+ (- v2 v1 -1) (* 2 halfsize) 1))))
  353.     (extract-canny-edges-and-grads 
  354.        subimage 
  355.        :mask-size mask-size
  356.        :lo-threshold lo-threshold
  357.        :hi-threshold hi-threshold)))))
  358.  
  359.  
  360. (defun find-oriented-2dworld-lines (edges orientation startx starty endx endy)
  361.   "  Edges are computed previously, for example by compute-2dworld-canny-edges.
  362.   Note: bounding box startx starty endx endy is given in 2d-world coordinates."
  363.   (when edges
  364.    (let ((image (edgeimg-image edges)))
  365.    (multiple-value-bind (u1 v1 u2 v2) 
  366.            (transform-2dworld-bounding-box image startx starty endx endy)
  367.      (setf orientation (transform-2dworld-angle image startx starty orientation))
  368.      (let* ((filtered (filter-edges-by-orientation edges (+ 90.0 orientation)
  369.                            :startx u1 :starty v1
  370.                            :endx u2 :endy v2))
  371.         (chains (extract-edgel-chains filtered 
  372.                       :startx u1 :starty v1
  373.                       :endx u2 :endy v2))
  374.         (lines (mapcar #'fit-line-to-pixlist chains)))
  375.        (transform-linelist-coordinates 
  376.     lines (image-to-2d-transform image)))))))
  377.  
  378. (defun put-lines-in-fs (linelist image-or-view fs)
  379.   (mapc #'(lambda (endpts)
  380.         (add-object 
  381.           (make-2d-curve
  382.             :vertices (make-vertex-array-from-vertex-list
  383.                 ;(list (list (car endpts) (cadr endpts) 0.0)
  384.                 ;      (list (third endpts) (fourth endpts) 0.0)))
  385.                 (list (list (car endpts) (cadr endpts) 0.0)
  386.                   (list (third endpts) (fourth endpts) 0.0)))
  387.         :closed-p t
  388.         :world (2d-world image-or-view))
  389.           fs))
  390.     linelist)
  391.   fs)
  392.  
  393.  
  394. (defun put-2.5d-lines-in-fs (linelist image-or-view fs)
  395.   (let ((trans (3d-to-2d-projection (2d-world image-or-view))))
  396.    (mapc #'(lambda (line)
  397.          (let ((x1y1 (epi::backproject-point
  398.                trans (car line) (cadr line) (fifth line)))
  399.            (x2y2 (epi::backproject-point 
  400.                trans (third line) (fourth line) (fifth line))))
  401.            (add-object 
  402.         (make-3d-curve
  403.          :vertices
  404.          (make-vertex-array-from-vertex-list
  405.           (list (list (car x1y1) (cadr x1y1) (fifth line))
  406.             (list (car x2y2) (cadr x2y2) (fifth line))))
  407.          :closed-p t
  408.          :world (3d-world image-or-view))
  409.         fs)))
  410.      linelist)
  411.    fs))
  412.  
  413.             
  414. #|
  415. (progn
  416.   (setf pane1 (ic::get-pane-named "pane-0-0"))
  417.   (setf pane2 (ic::get-pane-named "pane-1-0"))
  418.   (setf pane3 (ic::get-pane-named "pane-0-1"))
  419.   (setf pane4 (ic::get-pane-named "pane-1-1")))
  420.  
  421. (defun testlines (image fs orientation startx starty endx endy)
  422.   (put-lines-in-fs 
  423.      (find-oriented-lines image orientation
  424.     :startx startx :starty starty :endx endx :endy endy)
  425.      image fs))
  426. |#
  427.  
  428. #|
  429. (defun cle () 
  430.   (load (compile-file "/home/colossus/rcollins/topdown-line-finder")))
  431. |#
  432.