Skip to content

Commit 3036854

Browse files
Aravind-SureshAravind-Suresh
committed
Approach opencv#2 to Radon transform.
1 parent 8291657 commit 3036854

File tree

2 files changed

+156
-97
lines changed

2 files changed

+156
-97
lines changed

‎modules/ximgproc/include/opencv2/ximgproc.hpp‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -107,10 +107,8 @@ CV_EXPORTS_W void niBlackThreshold( InputArray _src, OutputArray _dst,
107107
* [radonTransform Computes the Radon Transform of the given input array]
108108
* @param src [Input image]
109109
* @param dst [Output sinogram, referenced]
110-
* @param accuracy [Minimum accuracy for angles ( IN DEGREES ), defaults to 1 degree]
111110
*/
112-
CV_EXPORTS_W void radonTransform(InputArray src, OutputArray dst,
113-
float accuracy = 1);
111+
CV_EXPORTS_W void radonTransform(InputArray src, OutputArray dst);
114112

115113
}
116114
}

‎modules/ximgproc/src/radon_transform.cpp‎

Lines changed: 155 additions & 94 deletions
Original file line numberDiff line numberDiff line change
@@ -65,126 +65,187 @@ namespace cv {
6565
#endif
6666

6767
/**
68-
* [linspace - Constructs linearly spaced vector of length n between start and end values ( both inclusive )]
69-
* @param start [start value]
70-
* @param end [end value]
71-
* @param n [length of the vector]
72-
* @return vec [linearly spaced vector]
68+
* Flags for range of angles
69+
* AR_x_y : theta lies in [x, y] ( both x and y inclusive )
7370
*/
74-
static inline std::vector<float> linspace(float start, float end, int n) {
75-
/**
76-
* A utility function for generating a vector of linearly spaced integers.
77-
*/
7871

79-
// Finding diff
80-
float diff = (end - start)*1.0f/(n-1);
81-
std::vector<float> vec(n);
82-
float val = start;
72+
enum ARFlags { AR_1_45, AR_46_90, AR_91_135, AR_136_179, AR_180 };
8373

84-
for(int i=0; i<n; ++i) {
85-
vec[i] = val;
86-
// Adding diff to val
87-
val += diff;
88-
}
74+
/**
75+
* [trans Function for computing radon transform given a particular angle theta]
76+
* @param src [Source image]
77+
* @param dst [Destination image]
78+
* @param m [Half the columns in src]
79+
* @param n [Half the rows in src]
80+
* @param rhoMax [Maximum value of rho ( restricted due to finite dimension of the image )]
81+
* @param theta [Given value of theta]
82+
* @param angleRange [One of the angle ranges]
83+
* @param maxTheta [Maximum value of theta ( in degrees ), defaults to 180 degrees]
84+
*/
85+
86+
void trans(Mat src, Mat& dst, int m, int n, int rhoMax, int theta, int angleRange, int maxTheta = 180) {
87+
int rc = rhoMax/2;
88+
89+
// Computing cos(theta) and sin(theta)
90+
float cosTh = static_cast<float>(std::cos(theta*CV_PI/180.0));
91+
float sinTh = static_cast<float>(std::sin(theta*CV_PI/180.0));
92+
93+
// Slope of the line = tan(90 - theta)
94+
float a = -cosTh/sinTh;
95+
96+
switch(angleRange) {
8997

90-
// Returning the vector
91-
return vec;
98+
// Range: [1, 45]
99+
case AR_1_45:
100+
for(int r = 1; r <= rhoMax; ++r) {
101+
float rho = r - rc;
102+
float b = rho/sinTh;
103+
104+
float yMax = std::min((int)(-a*m + b), n-1),
105+
yMin = std::max((int)(a*m + b), -n);
106+
107+
for(int y = yMin; y <= yMax; ++y) {
108+
float x = (y - b)/a;
109+
int xFloor = static_cast<int>(std::floor(x));
110+
float xUp = x - xFloor;
111+
float xLow = 1.0f - xUp;
112+
x = xFloor;
113+
x = static_cast<float>(std::max((int)(x), -m));
114+
x = static_cast<float>(std::min((int)(x), m-2));
115+
116+
dst.at<float>(rhoMax - r, maxTheta - theta - 1) += xLow*src.at<float>(y + n, x + m)
117+
+ xUp*src.at<float>(y + n, x + m + 1);
118+
}
119+
}
120+
break;
121+
122+
// Range: [46, 90]
123+
case AR_46_90:
124+
for(int r = 1; r <= rhoMax; ++r) {
125+
float rho = r - rc;
126+
float b = rho/sinTh;
127+
128+
float xMax = std::min((int)((-n - b)/a), m-1),
129+
xMin = std::max((int)((n - b)/a), -m);
130+
131+
for(int x = xMin; x <= xMax; ++x) {
132+
float y = a*x + b;
133+
int yFloor = static_cast<int>(std::floor(y));
134+
float yUp = y - yFloor;
135+
float yLow = 1.0f - yUp;
136+
y = yFloor;
137+
y = static_cast<float>(std::max((int)(y), -n));
138+
y = static_cast<float>(std::min((int)(y), n-2));
139+
140+
dst.at<float>(rhoMax - r, maxTheta - theta - 1) += yLow*src.at<float>(y + n, x + m)
141+
+ yUp*src.at<float>(y + n + 1, x + m);
142+
}
143+
}
144+
break;
145+
146+
// Range: [91, 135]
147+
case AR_91_135:
148+
for(int r = 1; r <= rhoMax; ++r) {
149+
float rho = r - rc;
150+
float b = rho/sinTh;
151+
152+
float xMax = std::min((int)((n - b)/a), m-1),
153+
xMin = std::max((int)((-n - b)/a), -m);
154+
155+
for(int x = xMin; x <= xMax; ++x) {
156+
float y = a*x + b;
157+
int yFloor = static_cast<int>(std::floor(y));
158+
float yUp = y - yFloor;
159+
float yLow = 1.0f - yUp;
160+
y = yFloor;
161+
y = static_cast<float>(std::max((int)(y), -n));
162+
y = static_cast<float>(std::min((int)(y), n-2));
163+
164+
dst.at<float>(rhoMax - r, maxTheta - theta - 1) += yLow*src.at<float>(y + n, x + m)
165+
+ yUp*src.at<float>(y + n + 1, x + m);
166+
}
167+
}
168+
break;
169+
170+
// Range: [136, 179]
171+
case AR_136_179:
172+
for(int r = 1; r <= rhoMax; ++r) {
173+
float rho = r - rc;
174+
float b = rho/sinTh;
175+
176+
float yMax = std::min((int)(a*m + b), n-1),
177+
yMin = std::max((int)(-a*m + b), -n);
178+
179+
for(int y = yMin; y <= yMax; ++y) {
180+
float x = (y - b)/a;
181+
int xFloor = static_cast<int>(std::floor(x));
182+
float xUp = x - xFloor;
183+
float xLow = 1.0f - xUp;
184+
x = xFloor;
185+
x = static_cast<float>(std::max((int)(x), -m));
186+
x = static_cast<float>(std::min((int)(x), m-2));
187+
188+
dst.at<float>(rhoMax - r, maxTheta - theta - 1) += xLow*src.at<float>(y + n, x + m)
189+
+ xUp*src.at<float>(y + n, x + m + 1);
190+
}
191+
}
192+
break;
193+
194+
// theta = 180 degrees | Integrating along a vertical line
195+
case AR_180:
196+
int rhoOffset = static_cast<int>((rhoMax - src.cols)/2);
197+
for(int x = 1; x <= src.cols; ++x) {
198+
int r = x + rhoOffset;
199+
r = rhoMax - r;
200+
for(int y = 0; y < src.rows; ++y) {
201+
dst.at<float>(r, 179) += src.at<float>(y, x);
202+
}
203+
}
204+
break;
205+
}
92206
}
93207

94208
/**
95-
* [meshgrid Returns a rectangular grid 2D space]
96-
* @param matX [Repeating vector along X direction]
97-
* @param matY [Repeating vector along Y direction]
98-
* @param X [Final meshgrid matrix X]
99-
* @param Y [Final meshgrid matrix Y]
209+
* [getAngleRange Returns flag AR_x_y for given angle theta]
210+
* @param theta [Given angle theta]
211+
* @return [flag AR_x_y]
100212
*/
101-
static inline void meshgrid(const Mat &matX, const Mat &matY, Mat &X, Mat &Y) {
102-
/**
103-
* A utitlity function to compute the meshgrid matrices for given input vectors.
104-
*/
105-
repeat(matX.reshape(1,1), static_cast<int>(matY.total()), 1, X);
106-
repeat(matY.reshape(1,1).t(), 1, static_cast<int>(matX.total()), Y);
213+
214+
int getAngleRange(int theta) {
215+
if(theta >= 1 && theta <= 45) return AR_1_45;
216+
else if(theta >=46 && theta <=90) return AR_46_90;
217+
else if(theta >= 91 && theta <= 135) return AR_91_135;
218+
else if(theta >= 136 && theta <= 179) return AR_136_179;
219+
else return AR_180;
107220
}
108221

109222
/**
110223
* [radonTransform Computes the Radon Transform of the given input array]
111224
* @param src [Input image]
112225
* @param dst [Output sinogram, referenced]
113-
* @param accuracy [Minimum accuracy for angles ( IN DEGREES ), defaults to 1 degree]
114226
*/
115-
void radonTransform(InputArray src, OutputArray dst, float accuracy) {
227+
228+
void radonTransform(InputArray src, OutputArray& dst) {
229+
116230
/**
117-
* Implementation of the Radon Transform
118-
*
119-
* Example:
120-
* ::
121-
* cv::Mat img, imgRadTrans;
122-
* :::: img holds some image ::::
123-
*
124-
* radonTransform(img, imgRadTrans, 0.5);
125-
* :::: imgRadTrans contains the radon transform ::::
231+
* Computes the Radon Transform for the given source image.
126232
*/
127-
// Asserting accuracy to be positive
128-
CV_Assert(accuracy > 0);
129233

234+
CV_Assert(src.type() == CV_8UC1);
130235
Mat srcMat_ = src.getMat(), srcMat;
131-
CV_Assert(srcMat_.type() == CV_8UC1);
132236

133237
srcMat_.convertTo(srcMat, CV_32F, 1.0/255);
134-
135238
int rows = srcMat.rows, cols = srcMat.cols;
136-
float diag = (float)std::sqrt(rows*rows + cols*cols);
137-
138-
Mat srcMatPadded;
139-
int padWidth = (int)((std::ceil(diag - cols) + 2)/2);
140-
int padHeight = (int)((std::ceil(diag - rows) + 2)/2);
141-
copyMakeBorder(srcMat, srcMatPadded, padWidth, padHeight, padWidth, padHeight, BORDER_CONSTANT, 0);
239+
int m = cols/2, n = rows/2;
142240

143-
int n = (int)(180/accuracy) + 1;
144-
std::vector<float> vecTheta = linspace(0, 180, n);
145-
int dstCols = static_cast<int>(vecTheta.size()), dstRows = srcMatPadded.rows;
241+
float rhoMax = static_cast<float>(std::ceil(std::sqrt(rows*rows + cols*cols)));
242+
int maxTheta = 180;
146243

147-
dst.create(Size(dstCols, dstRows), CV_32F);
244+
dst.create(Size(maxTheta, rhoMax + 1), CV_32F);
148245
Mat dstMat = dst.getMat();
149246

150-
int num = srcMatPadded.cols;
151-
std::vector<float> vecVal = linspace(-1, 1, num);
152-
153-
Mat meshX, meshY, meshXX, meshYY;
154-
meshgrid(Mat(vecVal), Mat(vecVal), meshX, meshY);
155-
156-
for(int i = 0; i < dstCols; ++i) {
157-
float theta = (90 - vecTheta[i])*static_cast<float>(CV_PI)/180.0f;
158-
float cosTh = std::cos(theta), sinTh = std::sin(theta);
159-
160-
Mat temp(srcMatPadded.size(), CV_32F);
161-
162-
// Rotating meshgrid entries by an angle theta
163-
meshXX = cosTh*meshX - sinTh*meshY;
164-
meshYY = sinTh*meshX + cosTh*meshY;
165-
166-
// Clamping values to [-1, 1]
167-
for(int l = 0; l < meshXX.cols; ++l) {
168-
for(int m = 0; m < meshXX.rows; ++m) {
169-
if(meshXX.at<float>(m, l) < -1) meshXX.at<float>(m, l) = -1;
170-
else if(meshXX.at<float>(m, l) > 1) meshXX.at<float>(m, l) = 1;
171-
172-
if(meshYY.at<float>(m, l) < -1) meshYY.at<float>(m, l) = -1;
173-
else if(meshYY.at<float>(m, l) > 1) meshYY.at<float>(m, l) = 1;
174-
}
175-
}
176-
177-
// Normalizing meshXX and meshYY to proper indices in the image
178-
normalize(meshXX, meshXX, 0, srcMatPadded.cols - 1, NORM_MINMAX, -1);
179-
normalize(meshYY, meshYY, 0, srcMatPadded.cols - 1, NORM_MINMAX, -1);
180-
181-
// Remapping, so as to interpolate for non-integer points
182-
remap(srcMatPadded, temp, meshXX, meshYY, INTER_LINEAR, BORDER_CONSTANT, Scalar(0));
183-
184-
// Computing columnwise sum and hence the a column of the sinogram
185-
Mat dstColReduced;
186-
reduce(temp, dstColReduced, 1, REDUCE_SUM, -1);
187-
dstColReduced.copyTo(dstMat.col(i));
247+
for(int i = 1; i <= 180; ++i) {
248+
trans(srcMat, dstMat, m, n, rhoMax, i, getAngleRange(i), maxTheta);
188249
}
189250
}
190251

0 commit comments

Comments
 (0)