Skip to content

Commit 78cf68c

Browse files
author
Cristy
committed
Rewrite rotating calipers algorithm
1 parent f327aa8 commit 78cf68c

File tree

1 file changed

+119
-181
lines changed

1 file changed

+119
-181
lines changed

‎MagickCore/attribute.c‎

Lines changed: 119 additions & 181 deletions
Original file line numberDiff line numberDiff line change
@@ -1085,96 +1085,72 @@ MagickExport size_t GetImageDepth(const Image *image,ExceptionInfo *exception)
10851085
%
10861086
*/
10871087

1088-
static double getDistance(const PointInfo *p,const PointInfo *q,
1089-
const PointInfo *v)
1088+
typedef struct _CaliperInfo
10901089
{
10911090
double
1092-
gamma;
1091+
area,
1092+
width,
1093+
height,
1094+
projection;
10931095

1094-
/*
1095-
Shortest distance between p to vector v through point q.
1096-
*/
1097-
gamma=v->y*PerceptibleReciprocal(v->x);
1098-
return(fabs(p->y-(q->y-gamma*q->x)-gamma*p->x)*
1099-
PerceptibleReciprocal(sqrt(gamma*gamma+1.0)));
1100-
}
1096+
ssize_t
1097+
p,
1098+
q,
1099+
v;
1100+
} CaliperInfo;
11011101

1102-
static PointInfo getIntersection(const PointInfo *a,const PointInfo *b,
1103-
const PointInfo *c,const PointInfo *d)
1102+
static double getDistance(PointInfo *p,PointInfo *q)
11041103
{
1105-
PointInfo
1106-
m,
1107-
p,
1108-
vertex;
1109-
1110-
/*
1111-
Intersection: b passing through a and d passing c.
1112-
*/
1113-
m.x=b->y*PerceptibleReciprocal(b->x);
1114-
p.x=a->y-m.x*a->x;
1115-
m.y=d->y*PerceptibleReciprocal(d->x);
1116-
p.y=c->y-m.y*c->x;
1117-
vertex.x=(p.y-p.x)*PerceptibleReciprocal(m.x-m.y);
1118-
vertex.y=m.x*vertex.x+p.x;
1119-
return(vertex);
1104+
return((p->x-q->x)*(p->x-q->x)+(p->y-q->y)*(p->y-q->y));
11201105
}
11211106

1122-
static double getRadians(const PointInfo *p,const PointInfo *q)
1107+
static double getProjection(PointInfo *p,PointInfo *q,PointInfo *v)
11231108
{
11241109
double
1125-
gamma;
1110+
distance;
11261111

11271112
/*
1128-
Angle in radians between vector p and q.
1113+
Projection of vector (x,y) - p into a line passing through p and q.
11291114
*/
1130-
gamma=sqrt(p->x*p->x+p->y*p->y)*sqrt(q->x*q->x+q->y*q->y);
1131-
return(acos((p->x*q->x+p->y*q->y)*PerceptibleReciprocal(gamma)));
1132-
}
1133-
1134-
static PointInfo *getVertex(PointInfo *vertices,const ssize_t n,
1135-
const size_t number_vertices)
1136-
{
1137-
return(vertices+(n % number_vertices));
1115+
distance=getDistance(p,q);
1116+
if (distance < MagickEpsilon)
1117+
return(INFINITY);
1118+
return((q->x-p->x)*(v->x-p->x)+(v->y-p->y)*(q->y-p->y))/sqrt(distance);
11381119
}
11391120

1140-
static PointInfo rotateVector(const PointInfo *p,const double radians)
1121+
static double getFeretDiameter(PointInfo *p,PointInfo *q,PointInfo *v)
11411122
{
1142-
PointInfo
1143-
vector;
1123+
double
1124+
distance;
11441125

11451126
/*
1146-
Rotates the vector p.
1127+
Distance from a point (x,y) to a line passing through p and q.
11471128
*/
1148-
vector.x=p->x*cos(radians)-p->y*sin(radians);
1149-
vector.y=p->x*sin(radians)+p->y*cos(radians);
1150-
return(vector);
1129+
distance=getDistance(p,q);
1130+
if (distance < MagickEpsilon)
1131+
return(INFINITY);
1132+
return((q->x-p->x)*(v->y-p->y)-(v->x-p->x)*(q->y-p->y))/sqrt(distance);
11511133
}
11521134

11531135
MagickExport PointInfo *GetImageMinimumBoundingBox(Image *image,
11541136
size_t *number_vertices,ExceptionInfo *exception)
11551137
{
1138+
CaliperInfo
1139+
caliper_info;
1140+
11561141
double
1157-
caliper_area = 0.0,
1158-
caliper_height = 0.0,
1159-
caliper_radians = 0.0,
1160-
caliper_width = 0.0,
1161-
radians = 0.0;
1142+
angle,
1143+
diameter;
11621144

11631145
PointInfo
11641146
*bounding_box,
1165-
caliper[4] = { { 1.0, 0.0 }, {-1.0, 0.0 }, { 0.0,-1.0 }, { 0.0, 1.0 } },
1166-
support[4][2] = { { { 0.0, 0.0 }, { 0.0, 0.0 } },
1167-
{ { 0.0, 0.0 }, { 0.0, 0.0 } },
1168-
{ { 0.0, 0.0 }, { 0.0, 0.0 } },
1169-
{ { 0.0, 0.0 }, { 0.0, 0.0 } } },
11701147
*vertices;
11711148

1172-
size_t
1173-
hull_vertices;
1149+
register ssize_t
1150+
i;
11741151

1175-
ssize_t
1176-
corner[4] = { 0, 0, 0, 0 },
1177-
n;
1152+
size_t
1153+
number_hull_vertices;
11781154

11791155
/*
11801156
Generate the minimum bounding box with the "Rotating Calipers" algorithm.
@@ -1184,7 +1160,7 @@ MagickExport PointInfo *GetImageMinimumBoundingBox(Image *image,
11841160
if (image->debug != MagickFalse)
11851161
(void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
11861162
*number_vertices=0;
1187-
vertices=GetImageConvexHull(image,&hull_vertices,exception);
1163+
vertices=GetImageConvexHull(image,&number_hull_vertices,exception);
11881164
if (vertices == (PointInfo *) NULL)
11891165
return((PointInfo *) NULL);
11901166
*number_vertices=4;
@@ -1195,135 +1171,97 @@ MagickExport PointInfo *GetImageMinimumBoundingBox(Image *image,
11951171
vertices=(PointInfo *) RelinquishMagickMemory(vertices);
11961172
return((PointInfo *) NULL);
11971173
}
1198-
for (n=1; n < (ssize_t) hull_vertices; n++)
1199-
{
1200-
if (vertices[n].y < vertices[corner[0]].y)
1201-
corner[0]=n;
1202-
if (vertices[n].y > vertices[corner[1]].y)
1203-
corner[1]=n;
1204-
if (vertices[n].x < vertices[corner[2]].x)
1205-
corner[2]=n;
1206-
if (vertices[n].x > vertices[corner[3]].x)
1207-
corner[3]=n;
1208-
}
1209-
while (radians < MagickPI)
1174+
caliper_info.area=2.0*image->columns*image->rows;
1175+
caliper_info.width=(double) image->columns+image->rows;
1176+
caliper_info.height=0.0;
1177+
caliper_info.projection=0.0;
1178+
caliper_info.p=(-1);
1179+
caliper_info.q=(-1);
1180+
caliper_info.v=(-1);
1181+
for (i=0; i < (ssize_t) number_hull_vertices; i++)
12101182
{
12111183
double
1212-
angle[4],
1213-
area,
1214-
height,
1215-
min_angle,
1216-
width;
1217-
1218-
PointInfo
1219-
edge[4];
1220-
1221-
edge[0].x=getVertex(vertices,corner[0]+1,hull_vertices)->x-
1222-
getVertex(vertices,corner[0],hull_vertices)->x;
1223-
edge[0].y=getVertex(vertices,corner[0]+1,hull_vertices)->y-
1224-
getVertex(vertices,corner[0],hull_vertices)->y;
1225-
edge[1].x=getVertex(vertices,corner[1]+1,hull_vertices)->x-
1226-
getVertex(vertices,corner[1],hull_vertices)->x;
1227-
edge[1].y=getVertex(vertices,corner[1]+1,hull_vertices)->y-
1228-
getVertex(vertices,corner[1],hull_vertices)->y;
1229-
edge[2].x=getVertex(vertices,corner[2]+1,hull_vertices)->x-
1230-
getVertex(vertices,corner[2],hull_vertices)->x;
1231-
edge[2].y=getVertex(vertices,corner[2]+1,hull_vertices)->y-
1232-
getVertex(vertices,corner[2],hull_vertices)->y;
1233-
edge[3].x=getVertex(vertices,corner[3]+1,hull_vertices)->x-
1234-
getVertex(vertices,corner[3],hull_vertices)->x;
1235-
edge[3].y=getVertex(vertices,corner[3]+1,hull_vertices)->y-
1236-
getVertex(vertices,corner[3],hull_vertices)->y;
1237-
angle[0]=getRadians(&edge[0],&caliper[0]);
1238-
angle[1]=getRadians(&edge[1],&caliper[1]);
1239-
angle[2]=getRadians(&edge[2],&caliper[2]);
1240-
angle[3]=getRadians(&edge[3],&caliper[3]);
1241-
if ((IsNaN(angle[0]) != MagickFalse) ||
1242-
(IsNaN(angle[1]) != MagickFalse) ||
1243-
(IsNaN(angle[2]) != MagickFalse) ||
1244-
(IsNaN(angle[3]) != MagickFalse))
1245-
break;
1246-
area=0.0;
1247-
min_angle=MagickMin(MagickMin(MagickMin(angle[0],angle[1]),angle[2]),
1248-
angle[3]);
1249-
caliper[0]=rotateVector(&caliper[0],min_angle);
1250-
caliper[1]=rotateVector(&caliper[1],min_angle);
1251-
caliper[2]=rotateVector(&caliper[2],min_angle);
1252-
caliper[3]=rotateVector(&caliper[3],min_angle);
1253-
if (fabs(angle[0]-min_angle) < MagickEpsilon)
1254-
{
1255-
width=getDistance(getVertex(vertices,corner[1],hull_vertices),
1256-
getVertex(vertices,corner[0],hull_vertices),&caliper[0]);
1257-
height=getDistance(getVertex(vertices,corner[3],hull_vertices),
1258-
getVertex(vertices,corner[2],hull_vertices),&caliper[2]);
1259-
}
1260-
else
1261-
if (fabs(angle[1]-min_angle) < MagickEpsilon)
1184+
area = 0.0,
1185+
max_projection = 0.0,
1186+
min_diameter = -1.0,
1187+
min_projection = 0.0;
1188+
1189+
register ssize_t
1190+
j,
1191+
k;
1192+
1193+
ssize_t
1194+
p = -1,
1195+
q = -1,
1196+
v = -1;
1197+
1198+
for (j=0; j < (ssize_t) number_hull_vertices; j++)
1199+
{
1200+
double
1201+
diameter;
1202+
1203+
diameter=fabs(getFeretDiameter(&vertices[i],
1204+
&vertices[(i+1) % number_hull_vertices],&vertices[j]));
1205+
if (min_diameter < diameter)
12621206
{
1263-
width=getDistance(getVertex(vertices,corner[0],hull_vertices),
1264-
getVertex(vertices,corner[1],hull_vertices),&caliper[1]);
1265-
height=getDistance(getVertex(vertices,corner[3],hull_vertices),
1266-
getVertex(vertices,corner[2],hull_vertices),&caliper[2]);
1267-
}
1268-
else
1269-
if (fabs(angle[2]-min_angle) < MagickEpsilon)
1270-
{
1271-
width=getDistance(getVertex(vertices,corner[1],hull_vertices),
1272-
getVertex(vertices,corner[0],hull_vertices),&caliper[0]);
1273-
height=getDistance(getVertex(vertices,corner[3],hull_vertices),
1274-
getVertex(vertices,corner[2],hull_vertices),&caliper[2]);
1275-
}
1276-
else
1277-
{
1278-
width=getDistance(getVertex(vertices,corner[1],hull_vertices),
1279-
getVertex(vertices,corner[0],hull_vertices),&caliper[0]);
1280-
height=getDistance(getVertex(vertices,corner[2],hull_vertices),
1281-
getVertex(vertices,corner[3],hull_vertices),&caliper[3]);
1282-
}
1283-
radians+=min_angle;
1284-
area=width*height;
1285-
if ((fabs(caliper_area) < MagickEpsilon) || (area <= caliper_area))
1207+
min_diameter=diameter;
1208+
p=i;
1209+
q=(i+1) % number_hull_vertices;
1210+
v=j;
1211+
}
1212+
}
1213+
for (k=0; k < (ssize_t) number_hull_vertices; k++)
1214+
{
1215+
double
1216+
projection;
1217+
1218+
/*
1219+
Rotating calipers.
1220+
*/
1221+
projection=getProjection(&vertices[p],&vertices[q],&vertices[k]);
1222+
min_projection=MagickMin(min_projection,projection);
1223+
max_projection=MagickMax(max_projection,projection);
1224+
}
1225+
area=min_diameter*(max_projection-min_projection);
1226+
if (caliper_info.area > area)
12861227
{
1287-
support[0][0]=(*getVertex(vertices,corner[0],hull_vertices));
1288-
support[0][1]=caliper[0];
1289-
support[1][0]=(*getVertex(vertices,corner[1],hull_vertices));
1290-
support[1][1]=caliper[1];
1291-
support[2][0]=(*getVertex(vertices,corner[2],hull_vertices));
1292-
support[2][1]=caliper[2];
1293-
support[3][0]=(*getVertex(vertices,corner[3],hull_vertices));
1294-
support[3][1]=caliper[3];
1295-
caliper_radians=radians;
1296-
caliper_area=area;
1297-
caliper_width=width;
1298-
caliper_height=height;
1299-
}
1300-
if (fabs(angle[0]-min_angle) < MagickEpsilon)
1301-
corner[0]++;
1302-
else
1303-
if (fabs(angle[1]-min_angle) < MagickEpsilon)
1304-
corner[1]++;
1305-
else
1306-
if (fabs(angle[2]-min_angle) < MagickEpsilon)
1307-
corner[2]++;
1308-
else
1309-
corner[3]++;
1228+
caliper_info.area=area;
1229+
caliper_info.width=min_diameter;
1230+
caliper_info.height=max_projection-min_projection;
1231+
caliper_info.projection=max_projection;
1232+
caliper_info.p=p;
1233+
caliper_info.q=q;
1234+
caliper_info.v=v;
1235+
}
13101236
}
1311-
bounding_box[0]=getIntersection(&support[0][0],&support[0][1],&support[3][0],
1312-
&support[3][1]);
1313-
bounding_box[1]=getIntersection(&support[3][0],&support[3][1],&support[1][0],
1314-
&support[1][1]);
1315-
bounding_box[2]=getIntersection(&support[1][0],&support[1][1],&support[2][0],
1316-
&support[2][1]);
1317-
bounding_box[3]=getIntersection(&support[2][0],&support[2][1],&support[0][0],
1318-
&support[0][1]);
1237+
/*
1238+
Initialize minimum bounding box.
1239+
*/
1240+
diameter=getFeretDiameter(&vertices[caliper_info.p],
1241+
&vertices[caliper_info.q],&vertices[caliper_info.v]);
1242+
angle=atan2(vertices[caliper_info.q].y-vertices[caliper_info.p].y,
1243+
vertices[caliper_info.q].x-vertices[caliper_info.p].x);
1244+
bounding_box[0].x=vertices[caliper_info.p].x+cos(angle)*
1245+
caliper_info.projection;
1246+
bounding_box[0].y=vertices[caliper_info.p].y+sin(angle)*
1247+
caliper_info.projection;
1248+
bounding_box[1].x=bounding_box[0].x+cos(angle+MagickPI/2.0)*diameter;
1249+
bounding_box[1].y=bounding_box[0].y+sin(angle+MagickPI/2.0)*diameter;
1250+
bounding_box[2].x=bounding_box[1].x+cos(angle)*(-caliper_info.height);
1251+
bounding_box[2].y=bounding_box[1].y+sin(angle)*(-caliper_info.height);
1252+
bounding_box[3].x=bounding_box[2].x+cos(angle+MagickPI/2.0)*(-diameter);
1253+
bounding_box[3].y=bounding_box[2].y+sin(angle+MagickPI/2.0)*(-diameter);
1254+
/*
1255+
Export minimum bounding box properties.
1256+
*/
13191257
(void) FormatImageProperty(image,"minimum-bounding-box:area","%.*g",
1320-
GetMagickPrecision(),caliper_area);
1258+
GetMagickPrecision(),caliper_info.area);
13211259
(void) FormatImageProperty(image,"minimum-bounding-box:width","%.*g",
1322-
GetMagickPrecision(),caliper_width);
1260+
GetMagickPrecision(),caliper_info.width);
13231261
(void) FormatImageProperty(image,"minimum-bounding-box:height","%.*g",
1324-
GetMagickPrecision(),caliper_height);
1262+
GetMagickPrecision(),caliper_info.height);
13251263
(void) FormatImageProperty(image,"minimum-bounding-box:angle","%.*g",
1326-
GetMagickPrecision(),90.0-RadiansToDegrees(caliper_radians));
1264+
GetMagickPrecision(),RadiansToDegrees(angle));
13271265
vertices=(PointInfo *) RelinquishMagickMemory(vertices);
13281266
return(bounding_box);
13291267
}

0 commit comments

Comments
 (0)