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