1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/PTdecode/CImg-1.3.0/examples/radon_transform.cpp Mon Aug 03 14:09:20 2009 +0100 1.3 @@ -0,0 +1,402 @@ 1.4 +/* 1.5 + # 1.6 + # File : radon_transform.cpp 1.7 + # ( C++ source file ) 1.8 + # 1.9 + # Description : An implementation of the Radon Transform. 1.10 + # This file is a part of the CImg Library project. 1.11 + # ( http://cimg.sourceforge.net ) 1.12 + # 1.13 + # Copyright : David G. Starkweather 1.14 + # ( starkdg@sourceforge.net - starkweatherd@cox.net ) 1.15 + # 1.16 + # License : CeCILL v2.0 1.17 + # ( http://www.cecill.info/licences/Licence_CeCILL_V2-en.html ) 1.18 + # 1.19 + # This software is governed by the CeCILL license under French law and 1.20 + # abiding by the rules of distribution of free software. You can use, 1.21 + # modify and/ or redistribute the software under the terms of the CeCILL 1.22 + # license as circulated by CEA, CNRS and INRIA at the following URL 1.23 + # "http://www.cecill.info". 1.24 + # 1.25 + # As a counterpart to the access to the source code and rights to copy, 1.26 + # modify and redistribute granted by the license, users are provided only 1.27 + # with a limited warranty and the software's author, the holder of the 1.28 + # economic rights, and the successive licensors have only limited 1.29 + # liability. 1.30 + # 1.31 + # In this respect, the user's attention is drawn to the risks associated 1.32 + # with loading, using, modifying and/or developing or reproducing the 1.33 + # software by the user in light of its specific status of free software, 1.34 + # that may mean that it is complicated to manipulate, and that also 1.35 + # therefore means that it is reserved for developers and experienced 1.36 + # professionals having in-depth computer knowledge. Users are therefore 1.37 + # encouraged to load and test the software's suitability as regards their 1.38 + # requirements in conditions enabling the security of their systems and/or 1.39 + # data to be ensured and, more generally, to use and operate it in the 1.40 + # same conditions as regards security. 1.41 + # 1.42 + # The fact that you are presently reading this means that you have had 1.43 + # knowledge of the CeCILL license and that you accept its terms. 1.44 + # 1.45 +*/ 1.46 + 1.47 +#define ROUNDING_FACTOR(x) (((x) >= 0) ? 0.5 : -0.5) 1.48 +#include <cmath> 1.49 +#include "CImg.h" 1.50 +using namespace cimg_library; 1.51 + 1.52 +// The lines below are necessary when using a non-standard compiler as visualcpp6. 1.53 +#ifdef cimg_use_visualcpp6 1.54 +#define std 1.55 +#endif 1.56 +#ifdef min 1.57 +#undef min 1.58 +#undef max 1.59 +#endif 1.60 + 1.61 +#ifndef cimg_imagepath 1.62 +#define cimg_imagepath "img/" 1.63 +#endif 1.64 + 1.65 +CImg<double> GaussianKernel(double rho); 1.66 +CImg<float> ApplyGaussian(CImg<unsigned char> im,double rho); 1.67 +CImg<unsigned char> RGBtoGrayScale(CImg<unsigned char> &im); 1.68 +int GetAngle(int dy,int dx); 1.69 +CImg<unsigned char> CannyEdges(CImg<float> im, double T1, double T2,bool doHysteresis); 1.70 +CImg<> RadonTransform(CImg<unsigned char> im,int N); 1.71 + 1.72 +int main(int argc,char **argv) { 1.73 + cimg_usage("Illustration of the Radon Transform"); 1.74 + 1.75 + const char *file = cimg_option("-f",cimg_imagepath "parrot_original.ppm","path and file name"); 1.76 + const double sigma = cimg_option("-r",1.0,"blur coefficient for gaussian low pass filter (lpf)"), 1.77 + thresh1 = cimg_option("-t1",0.50,"lower threshold for canny edge detector"), 1.78 + thresh2 = cimg_option("-t2",1.25,"upper threshold for canny edge detector");; 1.79 + const int N = cimg_option("-n",64,"number of angles to consider in the Radon transform - should be a power of 2"); 1.80 + 1.81 + //color to draw lines 1.82 + const unsigned char green[] = {0,255,0}; 1.83 + CImg<unsigned char> src(file); 1.84 + 1.85 + int rhomax = (int)std::sqrt((double)(src.dimx()*src.dimx() + src.dimy()*src.dimy()))/2; 1.86 + 1.87 + if (cimg::dialog(cimg::basename(argv[0]), 1.88 + "Instructions:\n" 1.89 + "Click on space bar or Enter key to display Radon transform of given image\n" 1.90 + "Click on anywhere in the transform window to display a \n" 1.91 + "corresponding green line in the original image\n", 1.92 + "Start", "Quit",0,0,0,0, 1.93 + src.get_resize(100,100,1,3),true)) std::exit(0); 1.94 + 1.95 + //retrieve a grayscale from the image 1.96 + CImg<unsigned char> grayScaleIm; 1.97 + if ((src.dimv() == 3) && (src.dimx() > 0) && (src.dimy() > 0) && (src.dimz() == 1)) 1.98 + grayScaleIm = (CImg<unsigned char>)src.get_pointwise_norm(0).quantize(255,false); 1.99 + else if ((src.dimv() == 1)&&(src.dimx() > 0) && (src.dimy() > 0) && (src.dimz() == 1)) 1.100 + grayScaleIm = src; 1.101 + else { // image in wrong format 1.102 + if (cimg::dialog(cimg::basename("wrong file format"), 1.103 + "Incorrect file format\n","OK",0,0,0,0,0, 1.104 + src.get_resize(100,100,1,3),true)) std::exit(0); 1.105 + } 1.106 + 1.107 + //blur the image with a Gaussian lpf to remove spurious edges (e.g. noise) 1.108 + CImg<float> blurredIm = ApplyGaussian(grayScaleIm,sigma); 1.109 + 1.110 + //use canny edge detection algorithm to get edge map of the image 1.111 + //- the threshold values are used to perform hysteresis in the edge detection process 1.112 + CImg<unsigned char> cannyEdgeMap = CannyEdges(blurredIm,thresh1,thresh2,false); 1.113 + CImg<unsigned char> radonImage = *(new CImg<unsigned char>(500,400,1,1,0)); 1.114 + 1.115 + //display the two windows 1.116 + CImgDisplay dispImage(src,"original image"); 1.117 + dispImage.move(CImgDisplay::screen_dimx()/8,CImgDisplay::screen_dimy()/8); 1.118 + CImgDisplay dispRadon(radonImage,"Radon Transform"); 1.119 + dispRadon.move(CImgDisplay::screen_dimx()/4,CImgDisplay::screen_dimy()/4); 1.120 + CImgDisplay dispCanny(cannyEdgeMap,"canny edges"); 1.121 + //start main display loop 1.122 + while (!dispImage.is_closed && !dispRadon.is_closed && 1.123 + !dispImage.is_keyQ && !dispRadon.is_keyQ && 1.124 + !dispImage.is_keyESC && !dispRadon.is_keyESC){ 1.125 + 1.126 + CImgDisplay::wait(dispImage,dispRadon); 1.127 + 1.128 + if (dispImage.is_keySPACE || dispRadon.is_keySPACE) 1.129 + { 1.130 + radonImage = (CImg<unsigned char>)RadonTransform(cannyEdgeMap,N).quantize(255,false).resize(500,400); 1.131 + radonImage.display(dispRadon); 1.132 + } 1.133 + 1.134 + //when clicking on dispRadon window, draw line in original image window 1.135 + if (dispRadon.button) 1.136 + { 1.137 + const double rho = dispRadon.mouse_y*rhomax/dispRadon.dimy(), 1.138 + theta = (dispRadon.mouse_x*N/dispRadon.dimx())*2*cimg::valuePI/N, 1.139 + x = src.dimx()/2 + rho*std::cos(theta), 1.140 + y = src.dimy()/2 + rho*std::sin(theta); 1.141 + const int x0 = (int)(x + 1000*std::cos(theta + cimg::valuePI/2)), 1.142 + y0 = (int)(y + 1000*std::sin(theta + cimg::valuePI/2)), 1.143 + x1 = (int)(x - 1000*std::cos(theta + cimg::valuePI/2)), 1.144 + y1 = (int)(y - 1000*std::sin(theta + cimg::valuePI/2)); 1.145 + src.draw_line(x0,y0,x1,y1,green,1.0f,0xF0F0F0F0).display(dispImage); 1.146 + } 1.147 + } 1.148 + return 0; 1.149 +} 1.150 +/** 1.151 + * PURPOSE: create a 5x5 gaussian kernel matrix 1.152 + * PARAM rho - gaussiam equation parameter (default = 1.0) 1.153 + * RETURN CImg<double> the gaussian kernel 1.154 + **/ 1.155 + 1.156 +CImg<double> GaussianKernel(double sigma = 1.0) 1.157 +{ 1.158 + CImg<double> resultIm(5,5,1,1,0); 1.159 + int midX = 3, midY = 3; 1.160 + cimg_forXY(resultIm,X,Y) 1.161 + { 1.162 + resultIm(X,Y) = std::ceil(256.0*(std::exp(-(midX*midX + midY*midY)/(2*sigma*sigma)))/(2*cimg::valuePI*sigma*sigma)); 1.163 + } 1.164 + return resultIm; 1.165 +} 1.166 +/* 1.167 + * PURPOSE: convolve a given image with the gaussian kernel 1.168 + * PARAM CImg<unsigned char> im - image to be convolved upon 1.169 + * PARAM double sigma - gaussian equation parameter 1.170 + * RETURN CImg<float> image resulting from the convolution 1.171 + * */ 1.172 +CImg<float> ApplyGaussian(CImg<unsigned char> im,double sigma) 1.173 +{ 1.174 + CImg<float> smoothIm(im.dimx(),im.dimy(),1,1,0); 1.175 + 1.176 + //make gaussian kernel 1.177 + CImg<float> gk = GaussianKernel(sigma); 1.178 + //apply gaussian 1.179 + 1.180 + CImg_5x5(I,int); 1.181 + cimg_for5x5(im,X,Y,0,0,I) 1.182 + { 1.183 + float sum = 0; 1.184 + sum += gk(0,0)*Ibb + gk(0,1)*Ibp + gk(0,2)*Ibc + gk(0,3)*Ibn + gk(0,4)*Iba; 1.185 + sum += gk(1,0)*Ipb + gk(1,1)*Ipp + gk(1,2)*Ipc + gk(1,3)*Ipn + gk(1,4)*Ipa; 1.186 + sum += gk(2,0)*Icb + gk(2,1)*Icp + gk(2,2)*Icc + gk(2,3)*Icn + gk(2,4)*Ica; 1.187 + sum += gk(3,0)*Inb + gk(3,1)*Inp + gk(3,2)*Inc + gk(3,3)*Inn + gk(3,4)*Ina; 1.188 + sum += gk(4,0)*Iab + gk(4,1)*Iap + gk(4,2)*Iac + gk(4,3)*Ian + gk(4,4)*Iaa; 1.189 + smoothIm(X,Y)= sum/256; 1.190 + } 1.191 + 1.192 + return smoothIm; 1.193 +} 1.194 +/** 1.195 + * PURPOSE: convert a given rgb image to a MxNX1 single vector grayscale image 1.196 + * PARAM: CImg<unsigned char> im - rgb image to convert 1.197 + * RETURN: CImg<unsigned char> grayscale image with MxNx1x1 dimensions 1.198 + **/ 1.199 + 1.200 +CImg<unsigned char> RGBtoGrayScale(CImg<unsigned char> &im) 1.201 +{ 1.202 + CImg<unsigned char> grayImage(im.dimx(),im.dimy(),im.dimz(),1,0); 1.203 + if (im.dimv() == 3) 1.204 + { cimg_forXYZ(im,X,Y,Z) 1.205 + { 1.206 + 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)); 1.207 + } 1.208 + } 1.209 + grayImage.quantize(255,false); 1.210 + return grayImage; 1.211 +} 1.212 +/** 1.213 + * PURPOSE: aux. function used by CannyEdges to quantize an angle theta given by gradients, dx and dy 1.214 + * into 0 - 7 1.215 + * PARAM: dx,dy - gradient magnitudes 1.216 + * RETURN int value between 0 and 7 1.217 + **/ 1.218 +int GetAngle(int dy,int dx) 1.219 +{ 1.220 + double angle = cimg::abs(std::atan2((double)dy,(double)dx)); 1.221 + if ((angle >= -cimg::valuePI/8)&&(angle <= cimg::valuePI/8))//-pi/8 to pi/8 => 0 1.222 + return 0; 1.223 + else if ((angle >= cimg::valuePI/8)&&(angle <= 3*cimg::valuePI/8))//pi/8 to 3pi/8 => pi/4 1.224 + return 1; 1.225 + else if ((angle > 3*cimg::valuePI/8)&&(angle <= 5*cimg::valuePI/8))//3pi/8 to 5pi/8 => pi/2 1.226 + return 2; 1.227 + else if ((angle > 5*cimg::valuePI/8)&&(angle <= 7*cimg::valuePI/8))//5pi/8 to 7pi/8 => 3pi/4 1.228 + return 3; 1.229 + 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 1.230 + return 4; 1.231 + else return 0; 1.232 +} 1.233 +/** 1.234 + * PURPOSE: create an edge map of the given image with hysteresis using thresholds T1 and T2 1.235 + * PARAMS: CImg<float> im the image to perform edge detection on 1.236 + * T1 lower threshold 1.237 + * T2 upper threshold 1.238 + * RETURN CImg<unsigned char> edge map 1.239 + **/ 1.240 +CImg<unsigned char> CannyEdges(CImg<float> im, double T1, double T2, bool doHysteresis=false) 1.241 +{ 1.242 + CImg<unsigned char> edges(im); 1.243 + CImg<float> secDerivs(im); 1.244 + secDerivs.fill(0); 1.245 + edges.fill(0); 1.246 + CImgList<float> gradients = im.get_gradient("xy",1); 1.247 + int image_width = im.dimx(); 1.248 + int image_height = im.dimy(); 1.249 + 1.250 + { cimg_forXY(im,X,Y) 1.251 + { 1.252 + double Gr = std::sqrt(std::pow((double)gradients[0](X,Y),2.0) + std::pow((double)gradients[1](X,Y),2.0)); 1.253 + double theta = GetAngle(Y,X); 1.254 + //if Gradient magnitude is positive and X,Y within the image 1.255 + //take the 2nd deriv in the appropriate direction 1.256 + if ((Gr > 0)&&(X < image_width-1)&&(Y < image_height - 1)) 1.257 + { 1.258 + if (theta == 0) 1.259 + secDerivs(X,Y) = im(X+2,Y) - 2*im(X+1,Y) + im(X,Y); 1.260 + else if (theta == 1) 1.261 + secDerivs(X,Y) = im(X+2,Y+2) - 2*im(X+1,Y+1) + im(X,Y); 1.262 + else if (theta == 2) 1.263 + secDerivs(X,Y) = im(X,Y+2) - 2*im(X,Y+1) + im(X,Y); 1.264 + else if (theta == 3) 1.265 + secDerivs(X,Y) = im(X+2,Y+2) - 2*im(X+1,Y+1) + im(X,Y); 1.266 + else if (theta == 4) 1.267 + secDerivs(X,Y) = im(X+2,Y) - 2*im(X+1,Y) + im(X,Y); 1.268 + } 1.269 + }} 1.270 + //for each 2nd deriv that crosses a zero point and magnitude passes the upper threshold. 1.271 + //Perform hysteresis in the direction of the gradient, rechecking the gradient 1.272 + //angle for each pixel that meets the threshold requirement. Stop checking when 1.273 + //the lower threshold is not reached. 1.274 + CImg_5x5(I,float); 1.275 + cimg_for5x5(secDerivs,X,Y,0,0,I) 1.276 + { 1.277 + if ( (Ipp*Ibb < 0)|| 1.278 + (Ipc*Ibc < 0)|| 1.279 + (Icp*Icb < 0) ) 1.280 + { 1.281 + double Gr = std::sqrt(std::pow((double)gradients[0](X,Y),2.0) + std::pow((double)gradients[1](X,Y),2.0)); 1.282 + int dir = GetAngle(Y,X); 1.283 + int Xt=X, Yt=Y, delta_x = 0, delta_y=0; 1.284 + double GRt = Gr; 1.285 + if (Gr >= T2) 1.286 + edges(X,Y) = 255; 1.287 + //work along the gradient in one direction 1.288 + if (doHysteresis) 1.289 + { 1.290 + while ((Xt > 0) && (Xt < image_width-1) && (Yt > 0) && (Yt < image_height-1)) 1.291 + { 1.292 + switch (dir){ 1.293 + case 0:delta_x=0;delta_y=1;break; 1.294 + case 1:delta_x=1;delta_y=1;break; 1.295 + case 2:delta_x=1;delta_y=0;break; 1.296 + case 3:delta_x=1;delta_y=-1;break; 1.297 + case 4:delta_x=0;delta_y=1;break; 1.298 + } 1.299 + Xt += delta_x; 1.300 + Yt += delta_y; 1.301 + GRt = std::sqrt(std::pow((double)gradients[0](Xt,Yt),2.0) + std::pow((double)gradients[1](Xt,Yt),2.0)); 1.302 + dir = GetAngle(Yt,Xt); 1.303 + if (GRt >= T1) 1.304 + edges(Xt,Yt) = 255; 1.305 + } 1.306 + //work along gradient in other direction 1.307 + Xt = X; Yt = Y; 1.308 + while ((Xt > 0) && (Xt < image_width-1) && (Yt > 0) && (Yt < image_height-1)) 1.309 + { 1.310 + switch (dir){ 1.311 + case 0:delta_x=0;delta_y=1;break; 1.312 + case 1:delta_x=1;delta_y=1;break; 1.313 + case 2:delta_x=1;delta_y=0;break; 1.314 + case 3:delta_x=1;delta_y=-1;break; 1.315 + case 4:delta_x=0;delta_y=1;break; 1.316 + } 1.317 + Xt -= delta_x; 1.318 + Yt -= delta_y; 1.319 + GRt = std::sqrt(std::pow((double)gradients[0](Xt,Yt),2.0) + std::pow((double)gradients[1](Xt,Yt),2.0)); 1.320 + dir = GetAngle(Yt,Xt); 1.321 + if (GRt >= T1) 1.322 + edges(Xt,Yt) = 255; 1.323 + } 1.324 + } 1.325 + } 1.326 + } 1.327 + return edges; 1.328 +} 1.329 +/** 1.330 + * PURPOSE: perform radon transform of given image 1.331 + * PARAM: CImg<unsigned char> im - image to detect lines 1.332 + * int N - number of angles to consider (should be a power of 2) 1.333 + * (the values of N will be spread over 0 to 2PI) 1.334 + * RETURN CImg<unsigned char> - transform of given image of size, N x D 1.335 + * D = rhomax = sqrt(dimx*dimx + dimy*dimy)/2 1.336 + **/ 1.337 +CImg<> RadonTransform(CImg<unsigned char> im,int N) 1.338 +{ 1.339 + int image_width = im.dimx(); 1.340 + int image_height = im.dimy(); 1.341 + 1.342 + //calc offsets to center the image 1.343 + float xofftemp = image_width/2.0f - 1; 1.344 + float yofftemp = image_height/2.0f - 1; 1.345 + int xoffset = (int)std::floor(xofftemp + ROUNDING_FACTOR(xofftemp)); 1.346 + int yoffset = (int)std::floor(yofftemp + ROUNDING_FACTOR(yofftemp)); 1.347 + float dtemp = (float)std::sqrt((double)(xoffset*xoffset + yoffset*yoffset)); 1.348 + int D = (int)std::floor(dtemp + ROUNDING_FACTOR(dtemp)); 1.349 + 1.350 + CImg<> imRadon(N,D,1,1,0); 1.351 + 1.352 + //for each angle k to consider 1.353 + for (int k= 0 ; k < N; k++) 1.354 + { 1.355 + //only consider from PI/8 to 3PI/8 and 5PI/8 to 7PI/8 1.356 + //to avoid computational complexity of a steep angle 1.357 + if (k == 0){k = N/8;continue;} 1.358 + else if (k == (3*N/8 + 1)){ k = 5*N/8;continue;} 1.359 + else if (k == 7*N/8 + 1){k = N; continue;} 1.360 + 1.361 + //for each rho length, determine linear equation and sum the line 1.362 + //sum is to sum the values along the line at angle k2pi/N 1.363 + //sum2 is to sum the values along the line at angle k2pi/N + N/4 1.364 + //The sum2 is performed merely by swapping the x,y axis as if the image were rotated 90 degrees. 1.365 + for (int d=0; d < D; d++){ 1.366 + double theta = 2*k*cimg::valuePI/N;//calculate actual theta 1.367 + double alpha = std::tan(theta+cimg::valuePI/2);//calculate the slope 1.368 + double beta_temp = -alpha*d*std::cos(theta) + d*std::sin(theta);//y-axis intercept for the line 1.369 + int beta = (int)std::floor(beta_temp + ROUNDING_FACTOR(beta_temp)); 1.370 + //for each value of m along x-axis, calculate y 1.371 + //if the x,y location is within the boundary for the respective image orientations, add to the sum 1.372 + unsigned int sum1 = 0, 1.373 + sum2 = 0; 1.374 + int M = (image_width >= image_height) ? image_width : image_height; 1.375 + for (int m=0;m < M; m++) 1.376 + { 1.377 + //interpolate in-between values using nearest-neighbor approximation 1.378 + //using m,n as x,y indices into image 1.379 + double n_temp = alpha*(m-xoffset) + beta; 1.380 + int n = (int)std::floor(n_temp + ROUNDING_FACTOR(n_temp)); 1.381 + if ((m < image_width) && (n + yoffset >= 0) && (n + yoffset < image_height)) 1.382 + { 1.383 + sum1 += im(m, n + yoffset); 1.384 + } 1.385 + n_temp = alpha*(m-yoffset) + beta; 1.386 + n = (int)std::floor(n_temp + ROUNDING_FACTOR(n_temp)); 1.387 + if ((m < image_height)&&(n + xoffset >= 0)&&(n + xoffset < image_width)) 1.388 + { 1.389 + sum2 += im(-(n + xoffset) + image_width - 1, m); 1.390 + } 1.391 + } 1.392 + //assign the sums into the result matrix 1.393 + imRadon(k,d) = (float)sum1; 1.394 + //assign sum2 to angle position for theta+PI/4 1.395 + imRadon(((k + N/4)%N),d) = (float)sum2; 1.396 + } 1.397 + } 1.398 + return imRadon; 1.399 +} 1.400 +/* references: 1.401 + * 1. See Peter Toft's thesis on the Radon transform: http://petertoft.dk/PhD/index.html 1.402 + * While I changed his basic algorithm, the main idea is still the same and provides an excellent explanation. 1.403 + * 1.404 + * */ 1.405 +