ExampleDrawConic.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/**
   @example ExampleDrawConic.cpp
   @relates Conic2D,  ImageDraw
   @brief  Tried to implement a new of drawing conics.
        Unfortunately it doesn't work for some cases:
        - hyperbolas
        - and all other conics that have two or more parts visible on screen
        :-(
   @ingroup g_examples
   @author MIP
*/

#include <iostream>
#include <sstream>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <Base/Image/Image.hh>
#include <Base/Image/ImageIO.hh>
#include <Base/ImageUtils/ImageDraw.hh>

#include <Base/Geometry/HomgPoint2D.hh>
#include <Base/Geometry/HomgPoint3D.hh>
#include <Geometry/RMatrix.hh>
#include <Geometry/Conic2D.hh>

using namespace BIAS;
using namespace std;



// takes a conic and two points and returns the point that is closer to the
// conic
// if distance for both points is equal, the first will be returned
inline HomgPoint2D ConicMinDist(const Conic2D &conic,
                                const HomgPoint2D &p1, const HomgPoint2D &p2) {
  if (conic.LocatePoint(p1) <= conic.LocatePoint(p2)) {
    return p1;
  }
  else {
    return p2;
  }
}


bool ConicGetStartPoint(const Conic2D &conic,
                        unsigned int imageWidth, unsigned int imageHeight,
                        HomgPoint2D &startPoint) {
  // since the conic equation is of an implicit form, we first need to find
  // a pixel that is on (or at least near to) the conic
  // this is done by comparing neighbouring pixel. e.g. if pixel (x, y) is
  // inside the conic and pixel (x+1, y) is outside, then the conic is
  // between these two pixels.
  // first we check the image borders, so we get conics that cross the image
  // border very fast. furthermore this necessary for the main loop, since the
  // iteration will run from an (arbitrary) start point until the start point
  // is reached again or the image border.

  // some vars are redundant here and could be re-used
  HomgPoint2D p(0, 0), p1(0, 0), p2(0, 0);
  p.Homogenize();
  p1.Homogenize();
  p2.Homogenize();
  register int initialSign;
  int initialSign1, initialSign2;
  double dist;

  // first check upper and lower borders of the image
  p1[1] = 0;
  p2[1] = imageHeight - 1;

  dist = conic.LocatePoint(p1);
  if (dist > 0.0)
    initialSign1 = 1;
  else if (dist < 0.0)
    initialSign1 = -1;
  else
    initialSign1 = 0;

  dist = conic.LocatePoint(p2);
  if (dist > 0.0)
    initialSign2 = 1;
  else if (dist < 0.0)
    initialSign2 = -1;
  else
    initialSign2 = 0;

  for (register unsigned int x = 1; x < imageWidth; x++) {
    // upper border
    p1[0] = x;
    dist = conic.LocatePoint(p1);
    if (dist < 0.0 && initialSign1 > 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, 0), p1);
      return true;
    }
    if (dist > 0.0 && initialSign1 < 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, 0), p1);
      return true;
    }
    if (dist == 0.0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, 0), p1);
      return true;
    }

    // lower border
    p2[0] = x;
    dist = conic.LocatePoint(p2);
    if (dist < 0.0 && initialSign2 > 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, imageHeight - 1), p2);
      return true;
    }
    if (dist > 0.0 && initialSign2 < 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, imageHeight - 1), p2);
      return true;
    }
    if (dist == 0.0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, imageHeight - 1), p2);
      return true;
    }
  }

  // first check left and right borders of the image
  p1[0] = 0;
  p2[0] = imageWidth - 1;

  dist = conic.LocatePoint(p1);
  if (dist > 0.0)
    initialSign1 = 1;
  else if (dist < 0.0)
    initialSign1 = -1;
  else
    initialSign1 = 0;

  dist = conic.LocatePoint(p2);
  if (dist > 0.0)
    initialSign2 = 1;
  else if (dist < 0.0)
    initialSign2 = -1;
  else
    initialSign2 = 0;

  for (register unsigned int y = 1; y < imageHeight; y++) {
    // upper border
    p1[1] = y;
    dist = conic.LocatePoint(p1);
    if (dist < 0.0 && initialSign1 > 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(0, y), p1);
      return true;
    }
    if (dist > 0.0 && initialSign1 < 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(0, y), p1);
      return true;
    }
    if (dist == 0.0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(0, y), p1);
      return true;
    }

    // lower border
    p2[1] = y;
    dist = conic.LocatePoint(p2);
    if (dist < 0.0 && initialSign2 > 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(imageWidth - 1, y), p2);
      return true;
    }
    if (dist > 0.0 && initialSign2 < 0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(imageWidth - 1, y), p2);
      return true;
    }
    if (dist == 0.0) {
      startPoint = ConicMinDist(conic, HomgPoint2D(imageWidth - 1, y), p2);
      return true;
    }
  }

  // if nothing found so far, check the whole image
  for (register unsigned int y = 0; y < imageHeight; y++) {
    p[0] = -1;
    p[1] = y;
    dist = conic.LocatePoint(p);
    if (dist > 0.0)
      initialSign = 1;
    else if (dist < 0.0)
      initialSign = -1;
    else
      initialSign = 0;

    for (register unsigned int x = 0; x < imageWidth; x++) {
      p[0] = x;
      p[1] = y;
      dist = conic.LocatePoint(p);
      if (dist < 0.0 && initialSign > 0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, y), p);
        return true;
      }
      if (dist > 0.0 && initialSign < 0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, y), p);
        return true;
      }
      if (dist == 0.0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x - 1, y), p);
        return true;
      }
    }
  }

  // now run in perpendicular direction to catch x-parallel lines that we missed
  // however, this may result in a lot of cach misses, and thus is slow
  for (register unsigned int x = 0; x < imageWidth; x++) {
    p[0] = x;
    p[1] = -1;
    dist = conic.LocatePoint(p);
    if (dist > 0.0)
      initialSign = 1;
    else if (dist < 0.0)
      initialSign = -1;
    else
      initialSign = 0 ;

    for (register unsigned int y = 0; y < imageWidth; y++) {
      p[0] = x;
      p[1] = y;
      dist = conic.LocatePoint(p);
      if (dist < 0.0 && initialSign > 0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x, y - 1), p);
        return true;
      }
      if (dist > 0.0 && initialSign < 0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x, y - 1), p);
        return true;
      }
      if (dist == 0.0) {
        startPoint = ConicMinDist(conic, HomgPoint2D(x, y - 1), p);
        return true;
      }
    }
  }

  return false;
}


void ConicShadePixel(Image<unsigned char>& image, const HomgPoint2D &p) {
  // prepare vars
  unsigned char *pImageData = image.GetImageData();

  const unsigned int x = (unsigned int)rint(p[0]);
  const unsigned int y = (unsigned int)rint(p[1]);

  unsigned int offset = y * image.GetWidth() * image.GetChannelCount() +
                        x * image.GetChannelCount();

  for (unsigned int c = 0; c < image.GetChannelCount(); c++) {
    pImageData[offset + c] = 255;
  }
}


bool ConicGetNextPoint(const Conic2D &conic,
                       const Image<unsigned char> &visitedPoints,
                       const HomgPoint2D &curPoint,
                       HomgPoint2D &nextPoint) {
  HomgPoint2D p(0, 0), q(0, 0);
  p.Homogenize();
  q.Homogenize();

  double distP, distQ;
  bool pUnvisited, qUnvisited;

  const unsigned char **pVisitedPoints = visitedPoints.GetImageDataArray();

  cout << "getting next point" << endl;

  // iterate over a 3x3 patch around curPoint
  for (int y = -1; y <= 1; y++) {
    for (int x = -1; x <= 1; x++) {
      p[0] = curPoint[0] + x;
      p[1] = curPoint[1] + y;

      // compare to right neighbour
      if (x < 1) {
        q[0] = p[0] + 1;
        q[1] = p[1];

        distP = conic.LocatePoint(p);
        distQ = conic.LocatePoint(q);

        pUnvisited =
          pVisitedPoints[(unsigned int)p[0]][(unsigned int)p[1]] == 0;
        qUnvisited =
          pVisitedPoints[(unsigned int)q[0]][(unsigned int)q[1]] == 0;

        if (distP == 0.0 && pUnvisited) {
          nextPoint = p;
          return true;
        }
        if (distQ == 0.0 && qUnvisited) {
          nextPoint = q;
          return true;
        }
        if ((distP < 0.0 && distQ > 0.0) || (distP > 0.0 && distQ < 0.0)) {
          if (pUnvisited && qUnvisited) {
            nextPoint = ConicMinDist(conic, p, q);
            return true;
          }
        }
      }

      // compare to lower neighbour
      if (y < 1) {
        q[0] = p[0];
        q[1] = p[1] + 1;

        distP = conic.LocatePoint(p);
        distQ = conic.LocatePoint(q);

        pUnvisited =
          pVisitedPoints[(unsigned int)p[0]][(unsigned int)p[1]] == 0;
        qUnvisited =
          pVisitedPoints[(unsigned int)q[0]][(unsigned int)q[1]] == 0;

        if (distP == 0.0 && pUnvisited) {
          nextPoint = p;
          return true;
        }
        if (distQ == 0.0 && qUnvisited) {
          nextPoint = q;
          return true;
        }
        if ((distP < 0.0 && distQ > 0.0) || (distP > 0.0 && distQ < 0.0)) {
          if (pUnvisited && qUnvisited) {
            nextPoint = ConicMinDist(conic, p, q);
            return true;
          }
        }
      }
    }
  }

  return false;
}


int DrawConic(Image<unsigned char>& image, Conic2D &conic) {
  HomgPoint2D lastPoint, curPoint;

  Image<unsigned char> visitedPoints(image.GetWidth(), image.GetHeight(), 1);
  visitedPoints.SetZero();
  unsigned char **pVisitedPoints = visitedPoints.GetImageDataArray();

  if (!ConicGetStartPoint(conic, image.GetWidth(), image.GetHeight(), curPoint)) {
    // conic is outside of image or - if it is an ellipse - the radius is
    // less than 1 pixel
    return 1;
  }

  cout << endl << "got start point " << curPoint << endl;

  ConicShadePixel(image, curPoint);
  pVisitedPoints[(unsigned int)curPoint[0]][(unsigned int)curPoint[1]] = 1;

  // iterate over conic
  lastPoint = curPoint;
  while (ConicGetNextPoint(conic, visitedPoints, lastPoint, curPoint)) {
    cout << "got next point " << curPoint << endl;

    ConicShadePixel(image, curPoint);
    pVisitedPoints[(unsigned int)curPoint[0]][(unsigned int)curPoint[1]] = 1;
    lastPoint = curPoint;
  }

  return 0;
}



/**
 * Simple example that draws some conics
 */
int main(int argc, char** argv) {
  Image<unsigned char> image(400, 400, 1);
  image.SetZero();

  Conic2D conic;
  conic.SetEllipse(HomgPoint2D(200, 200), M_PI, 50, 50);

  // kevin's conic plotter
  conic.Draw(image);
    ImageIO::Save("conicKK.mip", image);

  image.SetZero();

  // my conic plotter
  DrawConic(image, conic);
    ImageIO::Save("conicRW.mip", image);

  return 0;
}