1 // -*- c-basic-offset: 4 -*-
2 /** @file CleanCP.cpp
3 *
4 * @brief algorithms for remove control points by statistic method
5 *
6 * the algorithm is based on ptoclean by Bruno Postle
7 *
8 * @author Thomas Modes
9 *
10 * $Id$
11 *
12 */
13
14 /* This is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public
16 * License as published by the Free Software Foundation; either
17 * version 2 of the License, or (at your option) any later version.
18 *
19 * This software is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 * Lesser General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public
25 * License along with this software. If not, see
26 * <http://www.gnu.org/licenses/>.
27 *
28 */
29
30 #include "CleanCP.h"
31 #include <algorithms/optimizer/PTOptimizer.h>
32 #include "algorithms/basic/CalculateCPStatistics.h"
33 #include "hugin_base/panotools/PanoToolsUtils.h"
34
35 namespace HuginBase {
36
getCPoutsideLimit_pair(Panorama pano,AppBase::ProgressDisplay & progress,double n)37 UIntSet getCPoutsideLimit_pair(Panorama pano, AppBase::ProgressDisplay& progress, double n)
38 {
39 CPVector allCP=pano.getCtrlPoints();
40 unsigned int nrImg=pano.getNrOfImages();
41 PanoramaOptions opts=pano.getOptions();
42 //set projection to equrectangular for optimisation
43 opts.setProjection(PanoramaOptions::EQUIRECTANGULAR);
44 pano.setOptions(opts);
45 UIntSet CPtoRemove;
46
47 // do optimisation of all images pair
48 // after it remove cp with errors > median/mean + n*sigma
49 for (unsigned int image1=0; image1<nrImg-1; image1++)
50 {
51 const SrcPanoImage& img=pano.getImage(image1);
52 for (unsigned int image2=image1+1; image2<nrImg; image2++)
53 {
54 //do not check linked image pairs
55 if(img.YawisLinkedWith(pano.getImage(image2)))
56 continue;
57 UIntSet Images;
58 Images.clear();
59 Images.insert(image1);
60 Images.insert(image2);
61 Panorama clean=pano.getSubset(Images);
62 // pictures should contain at least 2 control points
63 if(clean.getNrOfCtrlPoints()>1)
64 {
65 // remove all horizontal and vertical control points
66 CPVector cpl = clean.getCtrlPoints();
67 CPVector newCP;
68 for (CPVector::const_iterator it = cpl.begin(); it != cpl.end(); ++it)
69 if (it->mode == ControlPoint::X_Y)
70 newCP.push_back(*it);
71 clean.setCtrlPoints(newCP);
72
73 // we need at least 3 cp to optimize 3 variables: yaw, pitch and roll
74 if(clean.getNrOfCtrlPoints()>3)
75 {
76 //optimize position and hfov
77 OptimizeVector optvec;
78 std::set<std::string> imgopt;
79 //imgopt.insert("v");
80 optvec.push_back(imgopt);
81 imgopt.insert("r");
82 imgopt.insert("p");
83 imgopt.insert("y");
84 optvec.push_back(imgopt);
85 clean.setOptimizeVector(optvec);
86 PTools::optimize(clean);
87 cpl.clear();
88 cpl=clean.getCtrlPoints();
89 //calculate statistic and determine limit
90 double min,max,mean,var;
91 CalculateCPStatisticsError::calcCtrlPntsErrorStats(clean,min,max,mean,var);
92 // if the standard deviation is bigger than the value, assume we have a lot of
93 // false cp, in this case take the mean value directly as limit
94 double limit = (sqrt(var) > mean) ? mean : (mean + n*sqrt(var));
95
96 //identify cp with big error
97 unsigned int index=0;
98 unsigned int cpcounter=0;
99 for (CPVector::const_iterator it = allCP.begin(); it != allCP.end(); ++it)
100 {
101 if(it->mode == ControlPoint::X_Y)
102 if((it->image1Nr==image1 && it->image2Nr==image2) ||
103 (it->image1Nr==image2 && it->image2Nr==image1))
104 {
105 if (cpl[index].error>limit)
106 CPtoRemove.insert(cpcounter);
107 index++;
108 };
109 cpcounter++;
110 };
111 };
112 };
113 if (!progress.updateDisplayValue())
114 {
115 return CPtoRemove;
116 };
117 };
118 };
119
120 return CPtoRemove;
121 };
122
getCPoutsideLimit(Panorama pano,double n,bool skipOptimisation,bool includeLineCp)123 UIntSet getCPoutsideLimit(Panorama pano, double n, bool skipOptimisation, bool includeLineCp)
124 {
125 UIntSet CPtoRemove;
126 if(skipOptimisation)
127 {
128 //calculate current cp errors
129 HuginBase::PTools::calcCtrlPointErrors(pano);
130 }
131 else
132 {
133 //optimize pano, after optimization pano contains the cp errors of the optimized project
134 SmartOptimise::smartOptimize(pano);
135 };
136 CPVector allCP=pano.getCtrlPoints();
137 if(!includeLineCp)
138 {
139 //remove all horizontal and vertical CP for calculation of mean and sigma
140 CPVector CPxy;
141 for (CPVector::const_iterator it = allCP.begin(); it != allCP.end(); ++it)
142 {
143 if(it->mode == ControlPoint::X_Y)
144 CPxy.push_back(*it);
145 };
146 pano.setCtrlPoints(CPxy);
147 };
148 //calculate mean and sigma
149 double min,max,mean,var;
150 CalculateCPStatisticsError::calcCtrlPntsErrorStats(pano,min,max,mean,var);
151 if(!includeLineCp)
152 {
153 pano.setCtrlPoints(allCP);
154 };
155 // if the standard deviation is bigger than the value, assume we have a lot of
156 // false cp, in this case take the mean value directly as limit
157 double limit = (sqrt(var) > mean) ? mean : (mean + n*sqrt(var));
158
159 //now determine all control points with error > limit
160 unsigned int index=0;
161 for (CPVector::const_iterator it = allCP.begin(); it != allCP.end(); ++it)
162 {
163 if(it->error > limit)
164 {
165 if(includeLineCp)
166 {
167 // all cp are treated the same
168 // no need for further checks
169 CPtoRemove.insert(index);
170 }
171 else
172 {
173 //check only normal cp
174 if(it->mode == ControlPoint::X_Y)
175 {
176 CPtoRemove.insert(index);
177 };
178 };
179 };
180 index++;
181 };
182
183 return CPtoRemove;
184 };
185
getCPinMasks(HuginBase::Panorama pano)186 UIntSet getCPinMasks(HuginBase::Panorama pano)
187 {
188 HuginBase::UIntSet cps;
189 HuginBase::CPVector cpList=pano.getCtrlPoints();
190 if(!cpList.empty())
191 {
192 for(unsigned int i=0;i<cpList.size();i++)
193 {
194 HuginBase::ControlPoint cp=cpList[i];
195 // ignore line control points
196 if(cp.mode!=HuginBase::ControlPoint::X_Y)
197 continue;
198 bool insideMask=false;
199 // check first image
200 // remark: we could also use pano.getImage(cp.image1Nr).isInside(vigra::Point2D(cp.x1,cp.y1))
201 // this would also check the crop rectangles/circles
202 // but it would require that the pano is correctly align, otherwise the positive masks
203 // would not correctly checked
204 HuginBase::MaskPolygonVector masks=pano.getImage(cp.image1Nr).getMasks();
205 if(!masks.empty())
206 {
207 unsigned int j=0;
208 while((!insideMask) && (j<masks.size()))
209 {
210 insideMask=masks[j].isInside(hugin_utils::FDiff2D(cp.x1,cp.y1));
211 j++;
212 };
213 };
214 // and now the second
215 if(!insideMask)
216 {
217 masks=pano.getImage(cp.image2Nr).getMasks();
218 if(!masks.empty())
219 {
220 unsigned int j=0;
221 while((!insideMask) && (j<masks.size()))
222 {
223 insideMask=masks[j].isInside(hugin_utils::FDiff2D(cp.x2,cp.y2));
224 j++;
225 };
226 };
227 };
228 if(insideMask)
229 cps.insert(i);
230 };
231 }
232 return cps;
233 };
234
235 } // namespace
236