1 //: \file
2 // \author Kevin de Souza
3 // \date 31 January 2008
4 // \brief Program to crop a 3D image down to a specified bounding box.
5
6 #include <exception>
7 #include <iostream>
8 #include <algorithm>
9 #ifdef _MSC_VER
10 # include "vcl_msvc_warnings.h"
11 #endif
12 #include "vul/vul_arg.h"
13 #include <mbl/mbl_log.h>
14 #include <vimt3d/vimt3d_add_all_loaders.h>
15 #include <vimt3d/vimt3d_load.h>
16 #include <vimt3d/vimt3d_save.h>
17 #include <vil3d/vil3d_save.h>
18 #include <vil3d/vil3d_new.h>
19 #include <vil3d/vil3d_crop.h>
20 #include "vul/vul_file.h"
21 #include "vul/vul_string.h"
22 #include "vgl/vgl_point_3d.h"
23 #include "vgl/vgl_box_3d.h"
24
25
26 //=========================================================================
27 // Static function to create a static logger when first required
28 //=========================================================================
logger()29 static mbl_logger& logger()
30 {
31 static mbl_logger l("mul.tools.crop_image_3d");
32 return l;
33 }
34
35
36 //========================================================================
37 // Actual main function
38 //========================================================================
main2(int argc,char * argv[])39 int main2(int argc, char*argv[])
40 {
41 unsigned i0=0, j0=0, k0=0;
42 unsigned ni=10, nj=10, nk=10;
43 double fx0=0, fy0=0, fz0=0;
44 double fx1=1, fy1=1, fz1=1;
45 double x0=0, y0=0, z0=0;
46 double x1=100, y1=100, z1=100;
47
48 vul_arg_base::set_help_description(
49 "DESCRIPTION:\n"
50 "Crop a 3D image down to a specified bounding box.\n"
51 "The specified bounding box may be expanded slightly to enclose whole voxels.\n"
52 );
53
54 // Parse the program arguments
55 vul_arg<std::string> img_src(nullptr, "input image filename");
56 vul_arg<std::string> img_dst(nullptr, "output image filename");
57 vul_arg<std::vector<unsigned> > bbi("-bbi", "bounding box in image coords (i0,j0,k0,i1,j1,k1)");
58 vul_arg<std::vector<double> > bbf("-bbf", "bounding box in image fraction e.g. 0.2,0.2,0.2,0.75,0.75,0.75");
59 vul_arg<std::vector<double> > bbw("-bbw", "bounding box in world coords (x0,y0,z0,x1,y1,z1)");
60 vul_arg<std::vector<double> > cw("-cw", "crop width in world distances (xlo_d,ylo_d,zlo_d,xhi_d,yhi_d,zhi_d)");
61 vul_arg<bool> use_mm("-mm", "World coords are in units of millimetres (default=metres)", false);
62 vul_arg<bool> ignore_bounds_errors("-ib", "adjust any bounding box values outside image to image edge", false);
63 vul_arg_parse(argc, argv);
64
65 // Log the program arguments
66 if (logger().level() >= mbl_logger::INFO)
67 {
68 MBL_LOG(INFO, logger(), "crop_image_3d: ");
69 MBL_LOG(INFO, logger(), " img_src: " << img_src());
70 MBL_LOG(INFO, logger(), " img_dst: " << img_dst());
71 MBL_LOG(INFO, logger(), " use_mm: " << use_mm());
72 if (bbi.set())
73 {
74 std::ostream& log = logger().log(mbl_logger::INFO); // construct a new log message
75 log << " bbi: "; bbi.print_value(log);
76 log << std::endl; // terminate log message
77 }
78 if (bbf.set())
79 {
80 std::ostream& log = logger().log(mbl_logger::INFO); // construct a new log message
81 log << " bbf: "; bbf.print_value(log);
82 log << std::endl; // terminate log message
83 }
84 if (bbw.set())
85 {
86 std::ostream& log = logger().log(mbl_logger::INFO); // construct a new log message
87 log << " bbw: "; bbw.print_value(log);
88 log << std::endl; // terminate log message
89 }
90 if (cw.set())
91 {
92 std::ostream& log = logger().log(mbl_logger::INFO); // construct a new log message
93 log << " cw: "; cw.print_value(log);
94 log << std::endl; // terminate log message
95 }
96 }
97
98 // Count the number of bbox options specified - should be exactly 1.
99 unsigned nbb=0;
100 if (bbi.set()) nbb++;
101 if (bbf.set()) nbb++;
102 if (bbw.set()) nbb++;
103 if (cw.set()) nbb++;
104 if (nbb!=1)
105 {
106 std::cerr << "ERROR: specify exactly 1 of the -bbi, -bbf, -bbw, or -cw options." << std::endl;
107 return 1;
108 }
109
110 // Validate the bbi argument - should be 6 unsigneds, specifying
111 // lower corner (i0,j0,k0) and upper corner (i1,j1,k1) of included voxels
112 if (bbi.set())
113 {
114 if (bbi().size() != 6)
115 {
116 std::cerr << "ERROR: -bbi argument should contain exactly 6 unsigneds\n";
117 return 1;
118 }
119
120 if (bbi()[0] >= bbi()[3] || bbi()[1] >= bbi()[4] || bbi()[2] >= bbi()[5])
121 {
122 std::cerr << "ERROR: -bbi argument should indicate the lower and upper corners of a 3D box with strictly positive width, height and depth\n";
123 return 1;
124 }
125
126 i0 = bbi()[0]; j0 = bbi()[1]; k0 = bbi()[2];
127 ni = bbi()[3]-i0+1; nj = bbi()[4]-j0+1; nk = bbi()[5]-k0+1;
128 }
129
130 // Validate the bbf argument - should be 6 doubles, specifying
131 // lower corner fraction (fx0,fy0,fz0) and upper corner (fx1,fy1,fz1) of included voxels
132 if (bbf.set())
133 {
134 if (bbf().size() != 6)
135 {
136 std::cerr << "ERROR: -bbf argument should contain exactly 6 floats\n";
137 return 1;
138 }
139
140 if (bbf()[0] >= bbf()[3] || bbf()[1] >= bbf()[4] || bbf()[2] >= bbf()[5])
141 {
142 std::cerr << "ERROR: -bbf argument should indicate the lower and upper corners of a 3D box with strictly positive width, height and depth within [0,1]\n";
143 return 1;
144 }
145 fx0 = bbf()[0]; fy0 = bbf()[1]; fz0 = bbf()[2];
146 fx1 = bbf()[3]; fy1 = bbf()[4]; fz1 = bbf()[5];
147 }
148
149 // Validate the bbw argument - should be 6 doubles, specifying
150 // lower corner (x0,y0,z0) and upper corner (x1,y1,z1) in world coords
151 if (bbw.set())
152 {
153 if (bbw().size() != 6)
154 {
155 std::cerr << "ERROR: -bbw argument should contain exactly 6 floats\n";
156 return 1;
157 }
158
159 if (bbw()[0] >= bbw()[3] || bbw()[1] >= bbw()[4] || bbw()[2] >= bbw()[5])
160 {
161 std::cerr << "ERROR: -bbw argument should indicate the lower and upper corners of a 3D box with strictly positive width, height and depth.\n";
162 return 1;
163 }
164 x0 = bbw()[0]; y0 = bbw()[1]; z0 = bbw()[2];
165 x1 = bbw()[3]; y1 = bbw()[4]; z1 = bbw()[5];
166 }
167 // Validate the cw argument - should be 6 doubles, specifying
168 // lower corner shift up (x0,y0,z0) and upper corner shift down (x1,y1,z1) in world distances
169 if (cw.set())
170 {
171 if (cw().size() != 6)
172 {
173 std::cerr << "ERROR: -cw argument should contain exactly 6 floats\n";
174 return 1;
175 }
176
177 if (cw()[0] < 0 || cw()[1] || cw()[2] < 0 || cw()[3]<0 || cw()[4]<0 || cw()[5]<0)
178 {
179 std::cerr << "ERROR: -cw argument should specify cropping widths on each face of the bounding box. Negative values are not allowed.\n";
180 return 1;
181 }
182
183 x0 = cw()[0]; y0 = cw()[1]; z0 = cw()[2];
184 x1 = cw()[3]; y1 = cw()[4]; z1 = cw()[5];
185 }
186
187 // Determine the output filetype
188 std::string filetype = vul_file::extension(img_dst());
189 vul_string_left_trim(filetype, ".");
190 if (filetype.empty()) filetype = "v3i";
191
192 vimt3d_add_all_loaders();
193
194 vil3d_image_resource_sptr ir = vil3d_load_image_resource(img_src().c_str());
195 if (!ir)
196 {
197 std::cerr << "ERROR: Failed to load input image resource\n";
198 return 1;
199 }
200 MBL_LOG(INFO, logger(), "Loaded input image_resource");
201
202 MBL_LOG(INFO, logger(), "Input image: ");
203 MBL_LOG(INFO, logger(), " ni: " << ir->ni());
204 MBL_LOG(INFO, logger(), " nj: " << ir->nj());
205 MBL_LOG(INFO, logger(), " nk: " << ir->nk());
206 MBL_LOG(INFO, logger(), " np: " << ir->nplanes());
207
208 vimt3d_transform_3d w2i = vimt3d_load_transform(ir, use_mm());
209 MBL_LOG(INFO, logger(), "Loaded input image transform");
210
211 if (bbf.set())
212 {
213 if (ignore_bounds_errors())
214 {
215 if (fx0 < 0) fx0=0;
216 if (fy0 < 0) fy0=0;
217 if (fz0 < 0) fz0=0;
218 if (fx0 > 1) fx1=1;
219 if (fy0 > 1) fy1=1;
220 if (fz0 > 1) fz1=1;
221 }
222
223 // Convert image fraction values to voxel numbers
224
225 // Round lower bounds down
226 i0 = static_cast<unsigned>(std::floor((ir->ni()-1)*fx0));
227 j0 = static_cast<unsigned>(std::floor((ir->nj()-1)*fy0));
228 k0 = static_cast<unsigned>(std::floor((ir->nk()-1)*fz0));
229 // Round upper bounds up
230 auto i1 = static_cast<unsigned>(std::ceil((ir->ni()-1)*fx1));
231 auto j1 = static_cast<unsigned>(std::ceil((ir->nj()-1)*fy1));
232 auto k1 = static_cast<unsigned>(std::ceil((ir->nk()-1)*fz1));
233 ni = i1 - i0 + 1;
234 nj = j1 - j0 + 1;
235 nk = k1 - k0 + 1;
236 }
237 if (bbw.set())
238 {
239 // Convert world coords values to voxel numbers
240 vgl_point_3d<double> imlo = w2i(vgl_point_3d<double>(x0,y0,z0));
241 vgl_point_3d<double> imhi = w2i(vgl_point_3d<double>(x1,y1,z1));
242 imhi.set(imhi.x()*0.999999, imhi.y()*0.999999, imhi.z()*0.999999);
243 if (ignore_bounds_errors())
244 {
245 imlo.set(std::max<double>(imlo.x(),0),
246 std::max<double>(imlo.y(),0),
247 std::max<double>(imlo.z(),0) );
248 imhi.set(std::min<double>(imhi.x(),ir->ni()-1),
249 std::min<double>(imhi.y(),ir->nj()-1),
250 std::min<double>(imhi.z(),ir->nk()-1) );
251 }
252 // Round lower bounds down
253 i0 = static_cast<unsigned>(std::floor(imlo.x()));
254 j0 = static_cast<unsigned>(std::floor(imlo.y()));
255 k0 = static_cast<unsigned>(std::floor(imlo.z()));
256 // Round upper bounds up
257 auto i1 = static_cast<unsigned>(std::ceil(imhi.x()));
258 auto j1 = static_cast<unsigned>(std::ceil(imhi.y()));
259 auto k1 = static_cast<unsigned>(std::ceil(imhi.z()));
260
261
262 ni = i1 - i0 + 1;
263 nj = j1 - j0 + 1;
264 nk = k1 - k0 + 1;
265 }
266 if (cw.set())
267 {
268
269 // Convert world coords values to voxel numbers
270 vgl_box_3d<double> bb;
271 bb.add(w2i.inverse()(0,0,0));
272 bb.add(w2i.inverse()(ir->ni(),ir->nj(),ir->nk()));
273 bb.set_min_x(bb.min_x()+x0);
274 bb.set_min_y(bb.min_y()+y0);
275 bb.set_min_z(bb.min_z()+z0);
276 bb.set_max_x(bb.max_x()-x1);
277 bb.set_max_y(bb.max_y()-y1);
278 bb.set_max_z(bb.max_z()-z1);
279
280 if (bb.is_empty())
281 {
282 std::cerr << "ERROR: -cw argument should specify cropping widths that do not overlap.\n";
283 return 1;
284 }
285
286
287 vgl_point_3d<double> imhi = w2i(bb.max_point());
288 imhi.set(imhi.x()*0.999999, imhi.y()*0.999999, imhi.z()*0.999999);
289 vgl_point_3d<double> imlo = w2i(bb.min_point());
290
291 // Round lower bounds down
292 i0 = static_cast<unsigned>(std::max(0.0,std::floor(imlo.x())));
293 j0 = static_cast<unsigned>(std::max(0.0,std::floor(imlo.y())));
294 k0 = static_cast<unsigned>(std::max(0.0,std::floor(imlo.z())));
295 // Round upper bounds up
296 auto i1 = static_cast<unsigned>(std::ceil(imhi.x()));
297 auto j1 = static_cast<unsigned>(std::ceil(imhi.y()));
298 auto k1 = static_cast<unsigned>(std::ceil(imhi.z()));
299
300
301 ni = i1 - i0;
302 nj = j1 - j0;
303 nk = k1 - k0;
304 }
305 if (ignore_bounds_errors())
306 {
307 if (i0+ni > ir->ni()) ni=ir->ni() - i0;
308 if (j0+nj > ir->nj()) nj=ir->nj() - j0;
309 if (k0+nk > ir->nk()) nk=ir->nk() - k0;
310 }
311
312
313 if (logger().level()>=mbl_logger::DEBUG)
314 {
315 vgl_point_3d<double> img_orig_wc = w2i.inverse()(vgl_point_3d<double>(0,0,0));
316 MBL_LOG(DEBUG, logger(), "Image origin in world coords: " << img_orig_wc);
317 }
318
319 if (i0 >= ir->ni() || j0 >= ir->nj() || k0 > ir->nk())
320 {
321 std::cerr << "ERROR: Crop region bbox lower corner is outside input image " <<
322 vul_file::strip_directory(img_src()) << "\n";
323 return 2;
324 }
325 if (i0+ni > ir->ni())
326 {
327 std::cerr << "WARNING: Crop region bbox upper corner i was outside input image " <<
328 vul_file::strip_directory(img_src()) << "; truncating to fit.\n";
329 ni = ir->ni()-i0;
330 }
331 if (j0+nj > ir->nj())
332 {
333 std::cerr << "WARNING: Crop region bbox upper corner j was outside input image " <<
334 vul_file::strip_directory(img_src()) << "; truncating to fit.\n";
335 nj = ir->nj()-j0;
336 }
337 if (k0+nk > ir->nk())
338 {
339 std::cerr << "WARNING: Crop region bbox upper corner k was outside input image " <<
340 vul_file::strip_directory(img_src()) << "; truncating to fit.\n";
341 nk = ir->nk()-k0;
342 }
343
344 // Get a view of the cropped region
345 vil3d_image_view_base_sptr ivbp = ir->get_copy_view(i0, ni, j0, nj, k0, nk);
346
347 // Create output image resource
348 vil3d_image_resource_sptr ir2 = vil3d_new_image_resource(
349 img_dst().c_str(), ni, nj, nk, ivbp->nplanes(), ivbp->pixel_format(), filetype.c_str());
350 if (!ir2)
351 {
352 std::cerr << "ERROR: Failed to create output image resource\n";
353 return 2;
354 }
355 MBL_LOG(INFO, logger(), "Created output image_resource");
356
357 // Modify transform to account for the implicit translation
358 vimt3d_transform_3d transl;
359 transl.set_translation(-double(i0), -double(j0), -double(k0));
360 w2i = transl*w2i;
361 vimt3d_save_transform(ir2, w2i, use_mm());
362
363 if (logger().level()>=mbl_logger::DEBUG)
364 {
365 vgl_point_3d<double> img_orig_wc = w2i.inverse()(vgl_point_3d<double>(0,0,0));
366 MBL_LOG(DEBUG, logger(), "Image origin in world coords: " << img_orig_wc);
367 }
368
369 // Write the image to file
370 bool succ = ir2->put_view(*ivbp);
371 if (!succ)
372 {
373 std::cerr << "ERROR: Failed to put_view into output image resource\n";
374 return 3;
375 }
376 MBL_LOG(INFO, logger(), "Copied cropped image to output image_resource");
377 MBL_LOG(INFO, logger(), "Output image: ");
378 MBL_LOG(INFO, logger(), " ni: " << ir2->ni());
379 MBL_LOG(INFO, logger(), " nj: " << ir2->nj());
380 MBL_LOG(INFO, logger(), " nk: " << ir2->nk());
381 MBL_LOG(INFO, logger(), " np: " << ir2->nplanes());
382
383 return 0;
384 }
385
386
387 //========================================================================
388 // Exception-handling wrapper around main function
389 //========================================================================
main(int argc,char * argv[])390 int main(int argc, char*argv[])
391 {
392 // Initialize the logger
393 mbl_logger::root().load_log_config_file();
394
395 try
396 {
397 main2(argc, argv);
398 }
399 catch (std::exception& e)
400 {
401 std::cout << "Caught exception " << e.what() << std::endl;
402 return 3;
403 }
404 catch (...)
405 {
406 std::cout << "Caught unknown exception " << std::endl;
407 return 3;
408 }
409
410 return 0;
411 }
412