libptouch / file revision

 /*

  #

  #  File        : radon_transform.cpp

  #                ( C++ source file )

  #

  #  Description : An implementation of the Radon Transform.

  #                This file is a part of the CImg Library project.

  #                ( http://cimg.sourceforge.net )

  #

  #  Copyright   : David G. Starkweather

  #                ( starkdg@sourceforge.net - starkweatherd@cox.net )

  #

  #  License     : CeCILL v2.0

  #                ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html )

  #

  #  This software is governed by the CeCILL  license under French law and

  #  abiding by the rules of distribution of free software.  You can  use,

  #  modify and/ or redistribute the software under the terms of the CeCILL

  #  license as circulated by CEA, CNRS and INRIA at the following URL

  #  "http://www.cecill.info".

  #

  #  As a counterpart to the access to the source code and  rights to copy,

  #  modify and redistribute granted by the license, users are provided only

  #  with a limited warranty  and the software's author,  the holder of the

  #  economic rights,  and the successive licensors  have only  limited

  #  liability.

  #

  #  In this respect, the user's attention is drawn to the risks associated

  #  with loading,  using,  modifying and/or developing or reproducing the

  #  software by the user in light of its specific status of free software,

  #  that may mean  that it is complicated to manipulate,  and  that  also

  #  therefore means  that it is reserved for developers  and  experienced

  #  professionals having in-depth computer knowledge. Users are therefore

  #  encouraged to load and test the software's suitability as regards their

  #  requirements in conditions enabling the security of their systems and/or

  #  data to be ensured and,  more generally, to use and operate it in the

  #  same conditions as regards security.

  #

  #  The fact that you are presently reading this means that you have had

  #  knowledge of the CeCILL license and that you accept its terms.

  #

 */

 #define ROUNDING_FACTOR(x) (((x) >= 0) ? 0.5 : -0.5)

 #include <cmath>

 #include "CImg.h"

 using namespace cimg_library;

 // The lines below are necessary when using a non-standard compiler as visualcpp6.

 #ifdef cimg_use_visualcpp6

 #define std

 #endif

 #ifdef min

 #undef min

 #undef max

 #endif

 #ifndef cimg_imagepath

 #define cimg_imagepath "img/"

 #endif

 CImg<double> GaussianKernel(double rho);

 CImg<float> ApplyGaussian(CImg<unsigned char> im,double rho);

 CImg<unsigned char> RGBtoGrayScale(CImg<unsigned char> &im);

 int GetAngle(int dy,int dx);

 CImg<unsigned char> CannyEdges(CImg<float> im, double T1, double T2,bool doHysteresis);

 CImg<> RadonTransform(CImg<unsigned char> im,int N);

 int main(int argc,char **argv) {

   cimg_usage("Illustration of the Radon Transform");

   const char *file = cimg_option("-f",cimg_imagepath "parrot_original.ppm","path and file name");

   const double sigma = cimg_option("-r",1.0,"blur coefficient for gaussian low pass filter (lpf)"),

     thresh1 = cimg_option("-t1",0.50,"lower threshold for canny edge detector"),

     thresh2 = cimg_option("-t2",1.25,"upper threshold for canny edge detector");;

   const int N = cimg_option("-n",64,"number of angles to consider in the Radon transform - should be a power of 2");

   //color to draw lines

   const unsigned char green[] = {0,255,0};

   CImg<unsigned char> src(file);

   int rhomax = (int)std::sqrt((double)(src.dimx()*src.dimx() + src.dimy()*src.dimy()))/2;

   if (cimg::dialog(cimg::basename(argv[0]),

                    "Instructions:\n"

                    "Click on space bar or Enter key to display Radon transform of given image\n"

                    "Click on anywhere in the transform window to display a \n"

                    "corresponding green line in the original image\n",

                    "Start", "Quit",0,0,0,0,

                    src.get_resize(100,100,1,3),true)) std::exit(0);

   //retrieve a grayscale from the image

   CImg<unsigned char> grayScaleIm;

   if ((src.dimv() == 3) && (src.dimx() > 0) && (src.dimy() > 0) && (src.dimz() == 1))

     grayScaleIm = (CImg<unsigned char>)src.get_pointwise_norm(0).quantize(255,false);

   else if ((src.dimv() == 1)&&(src.dimx() > 0) && (src.dimy() > 0) && (src.dimz() == 1))

     grayScaleIm = src;

   else { // image in wrong format

     if (cimg::dialog(cimg::basename("wrong file format"),

                      "Incorrect file format\n","OK",0,0,0,0,0,

                      src.get_resize(100,100,1,3),true)) std::exit(0);

   }

   //blur the image with a Gaussian lpf to remove spurious edges (e.g. noise)

   CImg<float> blurredIm = ApplyGaussian(grayScaleIm,sigma);

   //use canny edge detection algorithm to get edge map of the image

   //- the threshold values are used to perform hysteresis in the edge detection process

   CImg<unsigned char> cannyEdgeMap = CannyEdges(blurredIm,thresh1,thresh2,false);

   CImg<unsigned char> radonImage = *(new CImg<unsigned char>(500,400,1,1,0));

   //display the two windows

   CImgDisplay dispImage(src,"original image");

   dispImage.move(CImgDisplay::screen_dimx()/8,CImgDisplay::screen_dimy()/8);

   CImgDisplay dispRadon(radonImage,"Radon Transform");

   dispRadon.move(CImgDisplay::screen_dimx()/4,CImgDisplay::screen_dimy()/4);

   CImgDisplay dispCanny(cannyEdgeMap,"canny edges");

   //start main display loop

   while (!dispImage.is_closed && !dispRadon.is_closed &&

          !dispImage.is_keyQ && !dispRadon.is_keyQ &&

          !dispImage.is_keyESC && !dispRadon.is_keyESC){

     CImgDisplay::wait(dispImage,dispRadon);

     if (dispImage.is_keySPACE || dispRadon.is_keySPACE)

     {

       radonImage = (CImg<unsigned char>)RadonTransform(cannyEdgeMap,N).quantize(255,false).resize(500,400);

       radonImage.display(dispRadon);

     }

     //when clicking on dispRadon window, draw line in original image window

     if (dispRadon.button)

       {

         const double rho = dispRadon.mouse_y*rhomax/dispRadon.dimy(),

           theta = (dispRadon.mouse_x*N/dispRadon.dimx())*2*cimg::valuePI/N,

           x = src.dimx()/2 + rho*std::cos(theta),

           y = src.dimy()/2 + rho*std::sin(theta);

         const int x0 = (int)(x + 1000*std::cos(theta + cimg::valuePI/2)),

           y0 = (int)(y + 1000*std::sin(theta + cimg::valuePI/2)),

           x1 = (int)(x - 1000*std::cos(theta + cimg::valuePI/2)),

           y1 = (int)(y - 1000*std::sin(theta + cimg::valuePI/2));

         src.draw_line(x0,y0,x1,y1,green,1.0f,0xF0F0F0F0).display(dispImage);

       }

   }

   return 0;

 }

 /**

  * PURPOSE: create a 5x5 gaussian kernel matrix

  * PARAM rho - gaussiam equation parameter (default = 1.0)

  * RETURN CImg<double> the gaussian kernel

  **/

 CImg<double> GaussianKernel(double sigma = 1.0)

 {

   CImg<double> resultIm(5,5,1,1,0);

   int midX = 3, midY = 3;

   cimg_forXY(resultIm,X,Y)

     {

       resultIm(X,Y) = std::ceil(256.0*(std::exp(-(midX*midX + midY*midY)/(2*sigma*sigma)))/(2*cimg::valuePI*sigma*sigma));

     }

   return resultIm;

 }

 /*

  * PURPOSE: convolve a given image with the gaussian kernel

  * PARAM CImg<unsigned char> im - image to be convolved upon

  * PARAM double sigma - gaussian equation parameter

  * RETURN CImg<float> image resulting from the convolution

  * */

 CImg<float> ApplyGaussian(CImg<unsigned char> im,double sigma)

 {

   CImg<float> smoothIm(im.dimx(),im.dimy(),1,1,0);

   //make gaussian kernel

   CImg<float> gk = GaussianKernel(sigma);

   //apply gaussian

   CImg_5x5(I,int);

   cimg_for5x5(im,X,Y,0,0,I)

     {

       float sum = 0;

       sum += gk(0,0)*Ibb + gk(0,1)*Ibp + gk(0,2)*Ibc + gk(0,3)*Ibn + gk(0,4)*Iba;

       sum += gk(1,0)*Ipb + gk(1,1)*Ipp + gk(1,2)*Ipc + gk(1,3)*Ipn + gk(1,4)*Ipa;

       sum += gk(2,0)*Icb + gk(2,1)*Icp + gk(2,2)*Icc + gk(2,3)*Icn + gk(2,4)*Ica;

       sum += gk(3,0)*Inb + gk(3,1)*Inp + gk(3,2)*Inc + gk(3,3)*Inn + gk(3,4)*Ina;

       sum += gk(4,0)*Iab + gk(4,1)*Iap + gk(4,2)*Iac + gk(4,3)*Ian + gk(4,4)*Iaa;

       smoothIm(X,Y)= sum/256;

     }

   return smoothIm;

 }

 /**

  * PURPOSE: convert a given rgb image to a MxNX1 single vector grayscale image

  * PARAM: CImg<unsigned char> im - rgb image to convert

  * RETURN: CImg<unsigned char> grayscale image with MxNx1x1 dimensions

  **/

 CImg<unsigned char> RGBtoGrayScale(CImg<unsigned char> &im)

 {

   CImg<unsigned char> grayImage(im.dimx(),im.dimy(),im.dimz(),1,0);

   if (im.dimv() == 3)

     {   cimg_forXYZ(im,X,Y,Z)

       {

         grayImage(X,Y,Z,0) = (unsigned char)(0.299*im(X,Y,Z,0) + 0.587*im(X,Y,Z,1) + 0.114*im(X,Y,Z,2));

       }

     }

   grayImage.quantize(255,false);

   return grayImage;

 }

 /**

  * PURPOSE: aux. function used by CannyEdges to quantize an angle theta given by gradients, dx and dy

  *  into 0 - 7

  * PARAM: dx,dy - gradient magnitudes

  * RETURN int value between 0 and 7

  **/

 int GetAngle(int dy,int dx)

 {

   double angle = cimg::abs(std::atan2((double)dy,(double)dx));

   if ((angle >= -cimg::valuePI/8)&&(angle <= cimg::valuePI/8))//-pi/8 to pi/8 => 0

     return 0;

   else if ((angle >= cimg::valuePI/8)&&(angle <= 3*cimg::valuePI/8))//pi/8 to 3pi/8 => pi/4

     return 1;

   else if ((angle > 3*cimg::valuePI/8)&&(angle <= 5*cimg::valuePI/8))//3pi/8 to 5pi/8 => pi/2

     return 2;

   else if ((angle > 5*cimg::valuePI/8)&&(angle <= 7*cimg::valuePI/8))//5pi/8 to 7pi/8 => 3pi/4

     return 3;

   else if (((angle > 7*cimg::valuePI/8) && (angle <= cimg::valuePI))  || ((angle <= -7*cimg::valuePI/8)&&(angle >= -cimg::valuePI))) //-7pi/8 to -pi OR 7pi/8 to pi => pi

     return 4;

   else return 0;

 }

 /**

  * PURPOSE: create an edge map of the given image with hysteresis using thresholds T1 and T2

  * PARAMS: CImg<float> im the image to perform edge detection on

  *                 T1 lower threshold

  *         T2 upper threshold

  * RETURN CImg<unsigned char> edge map

  **/

 CImg<unsigned char> CannyEdges(CImg<float> im, double T1, double T2, bool doHysteresis=false)

 {

   CImg<unsigned char> edges(im);

   CImg<float> secDerivs(im);

   secDerivs.fill(0);

   edges.fill(0);

   CImgList<float> gradients = im.get_gradient("xy",1);

   int image_width = im.dimx();

   int image_height = im.dimy();

   { cimg_forXY(im,X,Y)

     {

       double Gr = std::sqrt(std::pow((double)gradients[0](X,Y),2.0) + std::pow((double)gradients[1](X,Y),2.0));

       double theta = GetAngle(Y,X);

       //if Gradient magnitude is positive and X,Y within the image

       //take the 2nd deriv in the appropriate direction

       if ((Gr > 0)&&(X < image_width-1)&&(Y < image_height - 1))

         {

           if (theta == 0)

             secDerivs(X,Y) = im(X+2,Y) - 2*im(X+1,Y) + im(X,Y);

           else if (theta == 1)

             secDerivs(X,Y) = im(X+2,Y+2) - 2*im(X+1,Y+1) + im(X,Y);

           else if (theta == 2)

             secDerivs(X,Y) = im(X,Y+2) - 2*im(X,Y+1) + im(X,Y);

           else if (theta == 3)

             secDerivs(X,Y) = im(X+2,Y+2) - 2*im(X+1,Y+1) + im(X,Y);

           else if (theta == 4)

             secDerivs(X,Y) = im(X+2,Y) - 2*im(X+1,Y) + im(X,Y);

         }

     }}

   //for each 2nd deriv that crosses a zero point and magnitude passes the upper threshold.

   //Perform hysteresis in the direction of the gradient, rechecking the gradient

   //angle for each pixel that meets the threshold requirement.  Stop checking when

   //the lower threshold is not reached.

   CImg_5x5(I,float);

   cimg_for5x5(secDerivs,X,Y,0,0,I)

     {

       if (   (Ipp*Ibb < 0)||

              (Ipc*Ibc < 0)||

              (Icp*Icb < 0)   )

         {

           double Gr = std::sqrt(std::pow((double)gradients[0](X,Y),2.0) + std::pow((double)gradients[1](X,Y),2.0));

           int dir = GetAngle(Y,X);

           int Xt=X, Yt=Y, delta_x = 0, delta_y=0;

           double GRt = Gr;

           if (Gr >= T2)

             edges(X,Y) = 255;

           //work along the gradient in one direction

           if (doHysteresis)

             {

               while ((Xt > 0) && (Xt < image_width-1) && (Yt > 0) && (Yt < image_height-1))

                 {

                   switch (dir){

                   case 0:delta_x=0;delta_y=1;break;

                   case 1:delta_x=1;delta_y=1;break;

                   case 2:delta_x=1;delta_y=0;break;

                   case 3:delta_x=1;delta_y=-1;break;

                   case 4:delta_x=0;delta_y=1;break;

                   }

                   Xt += delta_x;

                   Yt += delta_y;

                   GRt = std::sqrt(std::pow((double)gradients[0](Xt,Yt),2.0) + std::pow((double)gradients[1](Xt,Yt),2.0));

                   dir = GetAngle(Yt,Xt);

                   if (GRt >= T1)

                     edges(Xt,Yt) = 255;

                 }

               //work along gradient in other direction

               Xt = X; Yt = Y;

               while ((Xt > 0) && (Xt < image_width-1) && (Yt > 0) && (Yt < image_height-1))

                 {

                   switch (dir){

                   case 0:delta_x=0;delta_y=1;break;

                   case 1:delta_x=1;delta_y=1;break;

                   case 2:delta_x=1;delta_y=0;break;

                   case 3:delta_x=1;delta_y=-1;break;

                   case 4:delta_x=0;delta_y=1;break;

                   }

                   Xt -= delta_x;

                   Yt -= delta_y;

                   GRt = std::sqrt(std::pow((double)gradients[0](Xt,Yt),2.0) + std::pow((double)gradients[1](Xt,Yt),2.0));

                   dir = GetAngle(Yt,Xt);

                   if (GRt >= T1)

                     edges(Xt,Yt) = 255;

                 }

             }

         }

     }

   return edges;

 }

 /**

  * PURPOSE: perform radon transform of given image

  * PARAM: CImg<unsigned char> im - image to detect lines

  *                int N - number of angles to consider (should be a power of 2)

  *                (the values of N will be spread over 0 to 2PI)

  * RETURN CImg<unsigned char> - transform of given image of size, N x D

  *                              D = rhomax = sqrt(dimx*dimx + dimy*dimy)/2

  **/

 CImg<> RadonTransform(CImg<unsigned char> im,int N)

 {

   int image_width = im.dimx();

   int image_height = im.dimy();

   //calc offsets to center the image

   float xofftemp = image_width/2.0f - 1;

   float yofftemp = image_height/2.0f - 1;

   int xoffset = (int)std::floor(xofftemp + ROUNDING_FACTOR(xofftemp));

   int yoffset = (int)std::floor(yofftemp + ROUNDING_FACTOR(yofftemp));

   float dtemp = (float)std::sqrt((double)(xoffset*xoffset + yoffset*yoffset));

   int D = (int)std::floor(dtemp + ROUNDING_FACTOR(dtemp));

   CImg<> imRadon(N,D,1,1,0);

   //for each angle k to consider

   for (int k= 0 ; k < N; k++)

     {

       //only consider from PI/8 to 3PI/8 and 5PI/8 to 7PI/8

       //to avoid computational complexity of a steep angle

       if (k == 0){k = N/8;continue;}

       else if (k == (3*N/8 + 1)){       k = 5*N/8;continue;}

       else if (k == 7*N/8 + 1){k = N;   continue;}

       //for each rho length, determine linear equation and sum the line

       //sum is to sum the values along the line at angle k2pi/N

       //sum2 is to sum the values along the line at angle k2pi/N + N/4

       //The sum2 is performed merely by swapping the x,y axis as if the image were rotated 90 degrees.

       for (int d=0; d < D; d++){

         double theta = 2*k*cimg::valuePI/N;//calculate actual theta

         double alpha = std::tan(theta+cimg::valuePI/2);//calculate the slope

         double beta_temp = -alpha*d*std::cos(theta) + d*std::sin(theta);//y-axis intercept for the line

         int beta = (int)std::floor(beta_temp + ROUNDING_FACTOR(beta_temp));

         //for each value of m along x-axis, calculate y

         //if the x,y location is within the boundary for the respective image orientations, add to the sum

         unsigned int sum1 = 0,

           sum2 = 0;

         int M = (image_width >= image_height) ? image_width : image_height;

         for (int m=0;m < M; m++)

           {

             //interpolate in-between values using nearest-neighbor approximation

             //using m,n as x,y indices into image

             double n_temp = alpha*(m-xoffset) + beta;

             int n = (int)std::floor(n_temp + ROUNDING_FACTOR(n_temp));

             if ((m < image_width) && (n + yoffset >= 0) && (n + yoffset < image_height))

               {

                 sum1 += im(m, n + yoffset);

               }

             n_temp = alpha*(m-yoffset) + beta;

             n = (int)std::floor(n_temp + ROUNDING_FACTOR(n_temp));

             if ((m < image_height)&&(n + xoffset >= 0)&&(n + xoffset < image_width))

               {

                 sum2 += im(-(n + xoffset) + image_width - 1, m);

               }

           }

         //assign the sums into the result matrix

         imRadon(k,d) = (float)sum1;

         //assign sum2 to angle position for theta+PI/4

         imRadon(((k + N/4)%N),d) = (float)sum2;

       }

     }

   return imRadon;

 }

 /* references:

  * 1. See Peter Toft's thesis on the Radon transform: http://petertoft.dk/PhD/index.html

  * While I changed his basic algorithm, the main idea is still the same and provides an excellent explanation.

  *

  * */