00001
00002
00027 #include <hugin_config.h>
00028 #include <hugin_version.h>
00029 #include <fstream>
00030 #include <sstream>
00031
00032 #include <vigra/error.hxx>
00033 #include <vigra/impex.hxx>
00034 #include <exiv2/image.hpp>
00035 #include <exiv2/exif.hpp>
00036 #ifdef WIN32
00037 #include <getopt.h>
00038 #else
00039 #include <unistd.h>
00040 #endif
00041
00042 #include <appbase/ProgressDisplayOld.h>
00043 #include <nona/SpaceTransform.h>
00044 #include <photometric/ResponseTransform.h>
00045
00046 #include <hugin_basic.h>
00047
00048 #include <foreign/lensdb/PTLensDB.h>
00049
00050 #include <tiffio.h>
00051 #include <vigra_ext/MultiThreadOperations.h>
00052 #include <vigra_ext/ImageTransforms.h>
00053
00054
00055 using namespace std;
00056 using namespace vigra;
00057
00058 using namespace hugin_utils;
00059 using namespace HuginBase;
00060
00061
00062 template <class SrcImgType, class FlatImgType, class DestImgType>
00063 void correctImage(SrcImgType & srcImg,
00064 const FlatImgType & srcFlat,
00065 SrcPanoImage src,
00066 vigra_ext::Interpolator interpolator,
00067 DestImgType & destImg,
00068 bool doCrop,
00069 AppBase::MultiProgressDisplay & progress);
00070
00071 template <class PIXELTYPE>
00072 void correctRGB(SrcPanoImage & src, ImageImportInfo & info, const char * outfile,
00073 bool crop, const std::string & compression, AppBase::MultiProgressDisplay & progress);
00074
00075 bool getPTLensCoef(const char * fn, string cameraMaker, string cameraName,
00076 string lensName, float focalLength, vector<double> & coeff);
00077
00078
00079
00080
00081 static void usage(const char * name)
00082 {
00083 cerr << name << ": correct lens distortion, vignetting and chromatic abberation" << std::endl
00084 << "fulla version " << DISPLAY_VERSION << endl
00085 << std::endl
00086 << "Usage: " << name << " [options] inputfile(s) " << std::endl
00087 << " option are: " << std::endl
00088 << " -g a:b:c:d Radial distortion coefficient for all channels, (a, b, c, d)" << std::endl
00089 << " -b a:b:c:d Radial distortion coefficents for blue channel, (a, b, c, d)" << std::endl
00090 << " this is applied on top of the -g distortion coefficients," << endl
00091 << " use for TCA corr" << std::endl
00092 << " -r a:b:c:d Radial distortion coefficents for red channel, (a, b, c, d)" << std::endl
00093 << " this is applied on top of the -g distortion coefficients," << endl
00094 << " use for TCA corr" << std::endl
00095 << " -p Try to read radial distortion coefficients for green" << endl
00096 << " channel from PTLens database" << std::endl
00097 << " -m Canon Camera manufacturer, for PTLens database query" << std::endl
00098 << " -n Camera Camera name, for PTLens database query" << std::endl
00099 << " -l Lens Lens name, for PTLens database query" << std::endl
00100 << " if not specified, a list of possible lenses is displayed" << std::endl
00101 << " -d 50 specify focal length in mm, for PTLens database query" << std::endl
00102 << " -s do not rescale the image to avoid black borders." << std::endl
00103 << endl
00104 << " -f filename Vignetting correction by flatfield division" << std::endl
00105 << " I = I / c, c = flatfield / mean(flatfield)" << std::endl
00106 << " -c a:b:c:d radial vignetting correction by division:" << std::endl
00107 << " I = I / c, c = a + b*r^2 + c*r^4 + d*r^6" << std::endl
00108 << " -i value gamma of input data. used for gamma correction" << std::endl
00109 << " before and after flatfield correction" << std::endl
00110 << " -t n Number of threads that should be used during processing" << std::endl
00111 << " -h Display help (this text)" << std::endl
00112 << " -o name set output filename. If more than one image is given," << std::endl
00113 << " the name will be uses as suffix (default suffix: _corr)" << std::endl
00114 << " -e value Compression of the output files" << std::endl
00115 << " For jpeg output: 0-100" << std::endl
00116 << " For tiff output: PACKBITS, DEFLATE, LZW" << std::endl
00117 << " -x X:Y Horizontal and vertical shift" << std::endl
00118 << " -v Verbose" << std::endl;
00119 }
00120
00121
00122 int main(int argc, char *argv[])
00123 {
00124
00125 const char * optstring = "e:g:b:r:pm:n:l:d:sf:c:i:t:ho:x:v";
00126 int o;
00127
00128
00129 opterr = 0;
00130
00131 vector<double> vec4(4);
00132 bool doFlatfield = false;
00133 bool doVigRadial = false;
00134 bool doCropBorders = true;
00135 unsigned nThreads=1;
00136 unsigned verbose = 0;
00137
00138 std::string batchPostfix("_corr");
00139 std::string outputFile;
00140 std::string compression;
00141 bool doPTLens = false;
00142 std::string cameraMaker;
00143 std::string cameraName;
00144 std::string lensName;
00145 float focalLength=0;
00146 double gamma = 1.0;
00147 double shiftX = 0;
00148 double shiftY = 0;
00149
00150 SrcPanoImage c;
00151 while ((o = getopt (argc, argv, optstring)) != -1)
00152 switch (o) {
00153 case 'e':
00154 compression = optarg;
00155 break;
00156 case 'r':
00157 if (sscanf(optarg, "%lf:%lf:%lf:%lf", &vec4[0], &vec4[1], &vec4[2], &vec4[3]) != 4)
00158 {
00159 std::cerr << std::endl << "Error: invalid -r argument" << std::endl <<std::endl;
00160 usage(argv[0]);
00161 return 1;
00162 }
00163 c.setRadialDistortionRed(vec4);
00164
00165 break;
00166 case 'g':
00167 if (sscanf(optarg, "%lf:%lf:%lf:%lf", &vec4[0], &vec4[1], &vec4[2], &vec4[3]) != 4)
00168 {
00169 std::cerr << std::endl << "Error: invalid -g argument" << std::endl <<std::endl;
00170 usage(argv[0]);
00171 return 1;
00172 }
00173 c.setRadialDistortion(vec4);
00174
00175 break;
00176 case 'b':
00177 if (sscanf(optarg, "%lf:%lf:%lf:%lf", &vec4[0], &vec4[1], &vec4[2], &vec4[3]) != 4)
00178 {
00179 std::cerr << std::endl << "Error: invalid -b argument" << std::endl <<std::endl;
00180 usage(argv[0]);
00181 return 1;
00182 }
00183 c.setRadialDistortionBlue(vec4);
00184
00185 break;
00186 case 's':
00187 doCropBorders = false;
00188 break;
00189 case 'f':
00190 c.setFlatfieldFilename(optarg);
00191 doFlatfield = true;
00192 break;
00193 case 'i':
00194 gamma = atof(optarg);
00195 c.setGamma(gamma);
00196 break;
00197 case 'p':
00198 doPTLens = true;
00199 break;
00200 case 'm':
00201 cameraMaker = optarg;
00202 doPTLens = true;
00203 break;
00204 case 'n':
00205 cameraName = optarg;
00206 doPTLens = true;
00207 break;
00208 case 'l':
00209 lensName = optarg;
00210 doPTLens = true;
00211 break;
00212 case 'd':
00213 focalLength = atof(optarg);
00214 doPTLens = true;
00215 break;
00216 case 'c':
00217 if (sscanf(optarg, "%lf:%lf:%lf:%lf", &vec4[0], &vec4[1], &vec4[2], &vec4[3]) !=4)
00218 {
00219 std::cerr << std::endl << "Error: invalid -c argument" << std::endl <<std::endl;
00220 usage(argv[0]);
00221 return 1;
00222 }
00223 c.setRadialVigCorrCoeff(vec4);
00224 doVigRadial=true;
00225 break;
00226 case '?':
00227 case 'h':
00228 usage(argv[0]);
00229 return 0;
00230 case 't':
00231 nThreads = atoi(optarg);
00232 break;
00233 case 'o':
00234 outputFile = optarg;
00235 break;
00236 case 'x':
00237 if (sscanf(optarg, "%lf:%lf", &shiftX, &shiftY) != 2)
00238 {
00239 std::cerr << std::endl << "Error: invalid -x argument" << std::endl <<std::endl;
00240 usage(argv[0]);
00241 return 1;
00242 }
00243 c.setRadialDistortionCenterShift(FDiff2D(shiftX, shiftY));
00244 break;
00245 case 'v':
00246 verbose++;
00247 break;
00248 default:
00249 abort ();
00250 }
00251
00252 if (doVigRadial && doFlatfield) {
00253 std::cerr << std::endl << "Error: cannot use -f and -c at the same time" << std::endl <<std::endl;
00254 usage(argv[0]);
00255 return 1;
00256 }
00257
00258 SrcPanoImage::VignettingCorrMode vm=SrcPanoImage::VIGCORR_NONE;
00259
00260 if (doVigRadial)
00261 vm = SrcPanoImage::VIGCORR_RADIAL;
00262 if (doFlatfield)
00263 vm = SrcPanoImage::VIGCORR_FLATFIELD;
00264
00265 vm = (SrcPanoImage::VignettingCorrMode) (vm | SrcPanoImage::VIGCORR_DIV);
00266 c.setVigCorrMode(vm);
00267
00268 unsigned nFiles = argc - optind;
00269 if (nFiles == 0) {
00270 std::cerr << std::endl << "Error: No input file(s) specified" << std::endl <<std::endl;
00271 usage(argv[0]);
00272 return 1;
00273 }
00274
00275
00276 vector<string> inFiles;
00277 vector<string> outFiles;
00278 if (nFiles == 1) {
00279 if (outputFile.length() !=0) {
00280 inFiles.push_back(string(argv[optind]));
00281 outFiles.push_back(outputFile);
00282 } else {
00283 string name = string(argv[optind]);
00284 inFiles.push_back(name);
00285 string basen = stripExtension(name);
00286 outFiles.push_back(basen.append(batchPostfix.append(".").append(getExtension(name))));
00287 }
00288 } else {
00289
00290 if (outputFile.length() != 0) {
00291 batchPostfix = outputFile;
00292 }
00293 for (int i = optind; i < argc; i++) {
00294 string name = string(argv[i]);
00295 inFiles.push_back(name);
00296 outFiles.push_back(stripExtension(name) + batchPostfix + "." + getExtension(name));
00297 }
00298 }
00299
00300
00301
00302 TIFFSetWarningHandler(0);
00303
00304 AppBase::StreamMultiProgressDisplay pdisp(cout);
00305
00306 if (nThreads == 0) nThreads = 1;
00307 vigra_ext::ThreadManager::get().setNThreads(nThreads);
00308
00309 try {
00310 vector<string>::iterator outIt = outFiles.begin();
00311 for (vector<string>::iterator inIt = inFiles.begin(); inIt != inFiles.end() ; ++inIt, ++outIt)
00312 {
00313 if (verbose > 0) {
00314 cerr << "Correcting " << *inIt << " -> " << *outIt << endl;
00315 }
00316 c.setFilename(*inIt);
00317
00318
00319 vigra::ImageImportInfo info(inIt->c_str());
00320 const char * pixelType = info.getPixelType();
00321 int bands = info.numBands();
00322 int extraBands = info.numExtraBands();
00323
00324
00325 if (doPTLens) {
00326 if (getPTLensCoef(inIt->c_str(), cameraMaker.c_str(), cameraName.c_str(),
00327 lensName.c_str(), focalLength, vec4))
00328 {
00329 c.setRadialDistortion(vec4);
00330 } else {
00331 cerr << "Error: could not extract correction parameters from PTLens database" << endl;
00332 return 1;
00333 }
00334 }
00335 c.setSize(info.size());
00336
00337 if (bands == 3 || bands == 4 && extraBands == 1) {
00338
00339 if (strcmp(pixelType, "UINT8") == 0) {
00340 correctRGB<RGBValue<UInt8> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00341 }
00342 else if (strcmp(pixelType, "UINT16") == 0) {
00343 correctRGB<RGBValue<UInt16> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00344 } else if (strcmp(pixelType, "INT16") == 0) {
00345 correctRGB<RGBValue<Int16> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00346 } else if (strcmp(pixelType, "UINT32") == 0) {
00347 correctRGB<RGBValue<UInt32> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00348 } else if (strcmp(pixelType, "FLOAT") == 0) {
00349 correctRGB<RGBValue<float> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00350 } else if (strcmp(pixelType, "DOUBLE") == 0) {
00351 correctRGB<RGBValue<double> >(c, info, outIt->c_str(), doCropBorders, compression, pdisp);
00352 }
00353 } else {
00354 DEBUG_ERROR("unsupported depth, only 3 channel images are supported");
00355 throw std::runtime_error("unsupported depth, only 3 channels images are supported");
00356 return 1;
00357 }
00358 }
00359 } catch (std::exception & e) {
00360 cerr << "caught exception: " << e.what() << std::endl;
00361 return 1;
00362 }
00363 return 0;
00364 }
00365
00366
00367
00368
00374 template <class SrcImgType, class FlatImgType, class DestImgType>
00375 void correctImage(SrcImgType & srcImg,
00376 const FlatImgType & srcFlat,
00377 SrcPanoImage src,
00378 vigra_ext::Interpolator interpolator,
00379 DestImgType & destImg,
00380 bool doCrop,
00381 AppBase::MultiProgressDisplay & progress)
00382 {
00383 typedef typename SrcImgType::value_type SrcPixelType;
00384 typedef typename DestImgType::value_type DestPixelType;
00385
00386 typedef typename vigra::NumericTraits<SrcPixelType>::RealPromote RSrcPixelType;
00387
00388
00389 progress.pushTask(AppBase::ProgressTask("correcting image", ""));
00390
00391 vigra::Diff2D shiftXY(- roundi(src.getRadialDistortionCenterShift().x),
00392 - roundi(src.getRadialDistortionCenterShift().y));
00393
00394 if( (src.getVigCorrMode() & SrcPanoImage::VIGCORR_FLATFIELD)
00395 || (src.getVigCorrMode() & SrcPanoImage::VIGCORR_RADIAL) )
00396 {
00397 src.setResponseType(HuginBase::SrcPanoImage::RESPONSE_LINEAR);
00398 Photometric::InvResponseTransform<SrcPixelType,SrcPixelType> invResp(src);
00399 invResp.enforceMonotonicity();
00400 if (src.getVigCorrMode() & SrcPanoImage::VIGCORR_FLATFIELD) {
00401 invResp.setFlatfield(&srcFlat);
00402 }
00403 vigra_ext::transformImageSpatial(srcImageRange(srcImg), destImage(srcImg), invResp, vigra::Diff2D(0,0));
00404 }
00405
00406 double scaleFactor=1.0;
00407
00408 if (doCrop) {
00409 scaleFactor=Nona::estScaleFactorForFullFrame(src);
00410 DEBUG_DEBUG("Black border correction scale factor: " << scaleFactor);
00411 double sf=scaleFactor;
00412 vector<double> radGreen = src.getRadialDistortion();
00413 for (int i=0; i < 4; i++) {
00414 radGreen[3-i] *=sf;
00415 sf *=scaleFactor;
00416 }
00417 src.setRadialDistortion(radGreen);
00418 }
00419
00420
00421 BImage alpha(destImg.size());
00422
00423 vigra_ext::PassThroughFunctor<typename SrcPixelType::value_type> ptf;
00424
00425 if (src.getCorrectTCA())
00426 {
00427
00428
00429
00430
00431
00432
00433
00434 Nona::SpaceTransform transfr;
00435 transfr.InitRadialCorrect(src, 0);
00436 if (transfr.isIdentity()) {
00437 vigra::copyImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), RedAccessor<SrcPixelType>()),
00438 destIter(destImg.upperLeft(), RedAccessor<DestPixelType>()));
00439 } else {
00440 vigra_ext::transformImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), RedAccessor<SrcPixelType>()),
00441 destIterRange(destImg.upperLeft(), destImg.lowerRight(), RedAccessor<DestPixelType>()),
00442 destImage(alpha),
00443 shiftXY,
00444 transfr,
00445 ptf,
00446 false,
00447 vigra_ext::INTERP_SPLINE_16,
00448 progress);
00449 }
00450
00451 Nona::SpaceTransform transfg;
00452 transfg.InitRadialCorrect(src, 1);
00453 if (transfg.isIdentity()) {
00454 vigra::copyImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), GreenAccessor<SrcPixelType>()),
00455 destIter(destImg.upperLeft(), GreenAccessor<DestPixelType>()));
00456 } else {
00457 transformImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), GreenAccessor<SrcPixelType>()),
00458 destIterRange(destImg.upperLeft(), destImg.lowerRight(), GreenAccessor<DestPixelType>()),
00459 destImage(alpha),
00460 shiftXY,
00461 transfg,
00462 ptf,
00463 false,
00464 vigra_ext::INTERP_SPLINE_16,
00465 progress);
00466 }
00467
00468 Nona::SpaceTransform transfb;
00469 transfb.InitRadialCorrect(src, 2);
00470 if (transfb.isIdentity()) {
00471 vigra::copyImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), BlueAccessor<SrcPixelType>()),
00472 destIter(destImg.upperLeft(), BlueAccessor<DestPixelType>()));
00473 } else {
00474 transformImage(srcIterRange(srcImg.upperLeft(), srcImg.lowerRight(), BlueAccessor<SrcPixelType>()),
00475 destIterRange(destImg.upperLeft(), destImg.lowerRight(), BlueAccessor<DestPixelType>()),
00476 destImage(alpha),
00477 shiftXY,
00478 transfb,
00479 ptf,
00480 false,
00481 vigra_ext::INTERP_SPLINE_16,
00482 progress);
00483 }
00484 } else {
00485
00486 Nona::SpaceTransform transf;
00487 transf.InitRadialCorrect(src, 1);
00488 vector <double> radCoeff = src.getRadialDistortion();
00489 if (transf.isIdentity() ||
00490 (radCoeff[0] == 0.0 && radCoeff[1] == 0.0 && radCoeff[2] == 0.0 && radCoeff[3] == 1.0))
00491 {
00492 vigra::copyImage(srcImageRange(srcImg),
00493 destImage(destImg));
00494 } else {
00495 vigra_ext::PassThroughFunctor<SrcPixelType> ptfRGB;
00496 transformImage(srcImageRange(srcImg),
00497 destImageRange(destImg),
00498 destImage(alpha),
00499 shiftXY,
00500 transf,
00501 ptfRGB,
00502 false,
00503 vigra_ext::INTERP_SPLINE_16,
00504 progress);
00505 }
00506 }
00507 }
00508
00509
00510
00511 template <class PIXELTYPE>
00512 void correctRGB(SrcPanoImage & src, ImageImportInfo & info, const char * outfile,
00513 bool crop, const std::string & compression, AppBase::MultiProgressDisplay & progress)
00514 {
00515 vigra::BasicImage<RGBValue<float> > srcImg(info.size());
00516 vigra::BasicImage<PIXELTYPE> output(info.size());
00517 importImage(info, destImage(srcImg));
00518 FImage flatfield;
00519 if (src.getVigCorrMode() & SrcPanoImage::VIGCORR_FLATFIELD) {
00520 ImageImportInfo finfo(src.getFlatfieldFilename().c_str());
00521 flatfield.resize(finfo.size());
00522 importImage(finfo, destImage(flatfield));
00523 }
00524 correctImage(srcImg, flatfield, src, vigra_ext::INTERP_SPLINE_16, output, crop, progress);
00525 ImageExportInfo outInfo(outfile);
00526 outInfo.setICCProfile(info.getICCProfile());
00527 outInfo.setPixelType(info.getPixelType());
00528 if (compression.size() > 0) {
00529 outInfo.setCompression(compression.c_str());
00530 }
00531 exportImage(srcImageRange(output), outInfo);
00532 }
00533
00534
00535 bool getPTLensCoef(const char * fn, string cameraMaker, string cameraName,
00536 string lensName, float focalLength, vector<double> & coeff)
00537 {
00538 int verbose_flag = 1;
00539 const char * profilePath = getenv("PTLENS_PROFILE");
00540 if (profilePath == NULL)
00541 {
00542 cerr << "ERROR: " << endl
00543 << " You need to specify the location of \"profile.txt\"." << endl
00544 << " Please set the PTLENS_PROFILE environment variable, for example:" << endl
00545 << " PTLENS_PROFILE=$HOME/.ptlens/profile.txt" << endl;
00546 return false;
00547 }
00548
00549 PTLDB_DB * db = PTLDB_readDB(profilePath);
00550 if (! db) {
00551 fprintf(stderr,"Failed to read PTLens profile: %s\n", profilePath);
00552 return false;
00553 }
00554
00555
00556
00557 Exiv2::Image::AutoPtr image = Exiv2::ImageFactory::open(fn);
00558 assert (image.get() != 0);
00559 image->readMetadata();
00560
00561 Exiv2::ExifData &exifData = image->exifData();
00562 if (exifData.empty()) {
00563 std::cout << fn << ": no EXIF data found in file" << std::endl;
00564 } else {
00565 if (cameraMaker.size() == 0) {
00566 cameraMaker = exifData["Exif.Image.Make"].toString();
00567 }
00568 if (cameraName.size() == 0) {
00569 cameraName = exifData["Exif.Image.Model"].toString();
00570 }
00571 if (focalLength == 0.0f) {
00572 focalLength = exifData["Exif.Photo.FocalLength"].toFloat();
00573 }
00574 }
00575
00576
00577
00578 PTLDB_CamNode * thisCamera = PTLDB_findCamera(db, cameraMaker.c_str(), cameraName.c_str());
00579 if (!thisCamera) {
00580 fprintf(stderr, "could not find camera: %s, %s\n", cameraMaker.c_str(), cameraName.c_str());
00581 return false;
00582 }
00583 PTLDB_LnsNode * thisLens = PTLDB_findLens(db, lensName.c_str(), thisCamera);
00584 if (thisLens == NULL)
00585 {
00586 fprintf(stderr, "Lens \"%s\" not found in database.\n", lensName.c_str());
00587 fprintf(stderr,"Available lenses for camera: %s\n", thisCamera->menuModel);
00588 PTLDB_LnsNode * lenses = PTLDB_findLenses(db, thisCamera);
00589 while (lenses != NULL)
00590 {
00591 fprintf(stderr,"%s\n", lenses->menuLens);
00592 lenses = lenses->nextLns;
00593 }
00594 return false;
00595 }
00596
00597 ImageImportInfo info(fn);
00598
00599 PTLDB_ImageInfo img;
00600 img.camera = thisCamera;
00601 img.lens = thisLens;
00602 img.width = info.size().x;
00603 img.height = info.size().y;
00604 img.focalLength = focalLength;
00605 img.converterDetected = 0;
00606 img.resize=0;
00607
00608 PTLDB_RadCoef coef;
00609 PTLDB_getRadCoefs(db, &img, &coef);
00610 if (verbose_flag)
00611 {
00612 fprintf(stderr,"%s %s, Lens %s @ %f mm\n", thisCamera->menuMake, thisCamera->menuModel, lensName.c_str(), focalLength);
00613 fprintf(stderr, "PTLens coeff: a=%8.6lf b=%8.6lf c=%8.6lf d=%8.6lf\n", coef.a, coef.b, coef.c, coef.d);
00614 }
00615 coeff[0] = coef.a;
00616 coeff[1] = coef.b;
00617 coeff[2] = coef.c;
00618 coeff[3] = coef.d;
00619 return true;
00620 }
00621