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