From: <st...@us...> - 2009-09-13 17:02:53
|
Revision: 4403 http://hugin.svn.sourceforge.net/hugin/?rev=4403&view=rev Author: stativ Date: 2009-09-13 17:02:46 +0000 (Sun, 13 Sep 2009) Log Message: ----------- Moving forward to the support of gray scale input. Finally this is one which shouldn't crash. Modified Paths: -------------- hugin/branches/gsoc2009_deghosting/src/deghosting/CMakeLists.txt hugin/branches/gsoc2009_deghosting/src/deghosting/deghosting_mask.cpp hugin/branches/gsoc2009_deghosting/src/deghosting/khan.h hugin/branches/gsoc2009_deghosting/src/tools/CMakeLists.txt hugin/branches/gsoc2009_deghosting/src/tools/hugin_hdrmerge.cpp Removed Paths: ------------- hugin/branches/gsoc2009_deghosting/src/deghosting/khan.cpp Modified: hugin/branches/gsoc2009_deghosting/src/deghosting/CMakeLists.txt =================================================================== --- hugin/branches/gsoc2009_deghosting/src/deghosting/CMakeLists.txt 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/deghosting/CMakeLists.txt 2009-09-13 17:02:46 UTC (rev 4403) @@ -1,5 +1,5 @@ -add_executable(deghosting_mask deghosting_mask.cpp deghosting.cpp khan.cpp) +add_executable(deghosting_mask deghosting_mask.cpp deghosting.cpp) target_link_libraries(deghosting_mask ${common_libs} ${image_libs}) install(TARGETS deghosting_mask DESTINATION ${BINDIR}) Modified: hugin/branches/gsoc2009_deghosting/src/deghosting/deghosting_mask.cpp =================================================================== --- hugin/branches/gsoc2009_deghosting/src/deghosting/deghosting_mask.cpp 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/deghosting/deghosting_mask.cpp 2009-09-13 17:02:46 UTC (rev 4403) @@ -344,7 +344,9 @@ Deghosting* deghoster = NULL; try{ - Khan khanDeghoster(inputFiles, flags, debugFlags, iterations, sigma, verbosity); + // FIXME + // handle gray input correctly + Khan<RGBValue<float> > khanDeghoster(inputFiles, flags, debugFlags, iterations, sigma, verbosity); deghoster = &khanDeghoster; deghoster->setCameraResponse(response); Deleted: hugin/branches/gsoc2009_deghosting/src/deghosting/khan.cpp =================================================================== --- hugin/branches/gsoc2009_deghosting/src/deghosting/khan.cpp 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/deghosting/khan.cpp 2009-09-13 17:02:46 UTC (rev 4403) @@ -1,405 +0,0 @@ - -/** - * Implementation of Khan's deghosting algorithm - * Copyright (C) 2009 Lukáš Jirkovský <l.j...@gm...> - * - * This program is free software; you can redistribute it and/or - * modify it under the terms of the GNU General Public License - * as published by the Free Software Foundation; either version 2 - * of the License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - *Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#include "khan.h" -#include <vigra/impex.hxx> -// for importImageAlpha, it's not really necessary -#include <vigra_ext/impexalpha.hxx> - -/*// for EMoR -#include <photometric/ResponseTransform.h> -// neded for zeroNegative(), used in response transformation -#include <vigra_ext/ImageTransforms.h> -// panotools format of image, needed for response transformation -#include <panotools/PanoToolsInterface.h>*/ - -#ifdef ATAN_KH -// neded for atan and abs -#include <cmath> -#endif - -// for GammaFunctor -#include <vigra/transformimage.hxx> - -// for resampleImage -#include <vigra/resizeimage.hxx> - -// for snprintf -#include <cstdio> -using std::snprintf; - -using std::endl; -using std::cout; - -namespace deghosting { - Khan::Khan(std::vector<std::string>& newInputFiles, const uint16_t newFlags, const uint16_t newDebugFlags, - int newIterations, double newSigma, int newVerbosity) { - try { - Deghosting::loadImages(newInputFiles); - Deghosting::setFlags(newFlags); - Deghosting::setDebugFlags(newDebugFlags); - Deghosting::setIterationNum(newIterations); - Deghosting::setVerbosity(newVerbosity); - - // I don't know why, but sigma for HDR input have to approximately 10 times smaller - // FIXME: Maybe it would be better to use different sigma for different images in case both HDR and LDR are mixed - const char * fileType= ImageImportInfo(newInputFiles[0].c_str()).getFileType(); - if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { - setSigma(newSigma/10); - } else { - setSigma(newSigma); - } - - for (unsigned int i=0; i<5; i++) - Deghosting::response.push_back(0); - PIPOW = sigma*std::sqrt(2*PI); - denom = 1/PIPOW; - } catch (...) { - throw; - } - } - - Khan::Khan(std::vector<ImageImportInfo>& newInputFiles, const uint16_t newFlags, const uint16_t newDebugFlags, - int newIterations, double newSigma, int newVerbosity) { - try { - Deghosting::loadImages(newInputFiles); - Deghosting::setFlags(newFlags); - Deghosting::setDebugFlags(newDebugFlags); - Deghosting::setIterationNum(newIterations); - Deghosting::setVerbosity(newVerbosity); - - // I don't know why, but sigma for HDR input have to approximately 10 times smaller - // FIXME: Maybe it would be better to use different sigma for different images in case both HDR and LDR are mixed - const char * fileType= newInputFiles[0].getFileType(); - if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { - setSigma(newSigma/10); - } else { - setSigma(newSigma); - } - - for (unsigned int i=0; i<5; i++) - Deghosting::response.push_back(0); - PIPOW = sigma*std::sqrt(2*PI); - denom = 1/PIPOW; - } catch (...) { - throw; - } - } - - void Khan::setSigma(double newSigma) { - sigma = newSigma; - } - - float Khan::Kh(deghosting::AlgTinyVector< float, 3 > x) { - #ifdef ATAN_KH - // good choice for sigma for this function is around 600 - return std::atan(-(x*x)+sigma)/PI + 0.5; - #else - // good choice for sigma for this function is around 30 - return (std::exp(-(x*x)/(2*sigma*sigma)) * denom); - #endif - } - - /*void Khan::linearizeRGB(std::string inputFile,FRGBImage *pInputImg) { - HuginBase::SrcPanoImage panoImg(inputFile); - panoImg.setResponseType(HuginBase::SrcPanoImage::RESPONSE_EMOR); - panoImg.setEMoRParams(response); - // response transform functor - HuginBase::Photometric::InvResponseTransform<RGBValue<float>, - RGBValue<float> > invResponse(panoImg); - invResponse.enforceMonotonicity(); - - // iterator to the upper left corner - FRGBImage::traverser imgIterSourceY = pInputImg->upperLeft(); - // iterator to he lower right corner - FRGBImage::traverser imgIterEnd = pInputImg->lowerRight(); - // iterator to the output - FRGBImage::traverser imgIterOut = pInputImg->upperLeft(); - // loop through the image - for (int y=1; imgIterSourceY.y != imgIterEnd.y; ++imgIterSourceY.y, ++imgIterOut.y, ++y) { - // iterator to the input - FRGBImage::traverser sx = imgIterSourceY; - // iterator to the output - FRGBImage::traverser dx = imgIterOut; - for (int x=1; sx.x != imgIterEnd.x; ++sx.x, ++dx.x, ++x) { - // transform image using response - *dx = vigra_ext::zeroNegative(invResponse(*sx, hugin_utils::FDiff2D(x, y))); - } - } - }*/ - - template <class PixelType> - void Khan::preprocessImage(unsigned int i, FImagePtr &weight, FLabImagePtr &LabImage) { - ImageImportInfo imgInfo(inputFiles[i]); - BasicImage<PixelType> * pInputImg = new BasicImage<PixelType>(imgInfo.size()); - BImage imgAlpha(imgInfo.size()); - weight = FImagePtr(new FImage(imgInfo.size())); - - // import alpha - importImageAlpha(imgInfo, destImage(*pInputImg),destImage(imgAlpha)); - - // linearize RGB and convert it to L*a*b image - //linearizeRGB(inputFiles[i], pInputImg); - - // take logarithm or gamma correction if the input images are HDR - // I'm not sure if it's the right way how to - // find out if they are HDR - const char * fileType= imgInfo.getFileType(); - if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { - // use gamma 2.2 - if (flags & ADV_GAMMA) { - // GammaFunctor is only in vigra 1.6 GRRR - // I have to use BrightnessContrastFunctor - // TODO: change to the GammaFunctor in the future - vigra::FindMinMax<float> minmax; - vigra::inspectImage(srcImageRange(*pInputImg), minmax); - transformImage(srcImageRange(*pInputImg),destImage(*pInputImg),BrightnessContrastFunctor<PixelType>(0.45,1.0,minmax.min, minmax.max)); - } else { - // take logarithm - transformImage(srcImageRange(*pInputImg),destImage(*pInputImg),LogarithmFunctor<PixelType>(1.0)); - } - } - - // generate initial weights using hat function - if (flags & ADV_NOINITWEIGHTS) { - // convert image to grayscale (and take logarithm) - RGBToGrayAccessor<FRGBImage::PixelType> color2gray; - transformImage(srcImageRange(*pInputImg, color2gray), destImage(*weight), log(Arg1()+Param(1.0f))); - } - else { - transformImage(srcImageRange(*pInputImg),destImage(*weight),HatFunctor<PixelType>()); - } - - // convert from linear RGB to L*a*b // - RGB2LabFunctor<float> RGB2Lab; - LabImage = FLabImagePtr(new FLabImage(imgInfo.size())); - transformImage(srcImageRange(*pInputImg), destImage(*LabImage), RGB2Lab); - - delete pInputImg; - pInputImg = 0; - } - - std::vector<FImagePtr> Khan::createWeightMasks() { - for (unsigned int i = 0; i < inputFiles.size(); i++) { - FImagePtr weight; - FLabImagePtr LabImage; - preprocessImage<RGBValue<float> >(i, weight, LabImage); - LabImages.push_back(LabImage); - weights.push_back(weight); - - // save init weights - if (debugFlags & SAVE_INITWEIGHTS) { - char tmpfn[100]; - snprintf(tmpfn, 99, "init_weights_%d.tiff", i); - ImageExportInfo exWeights(tmpfn); - exportImage(srcImageRange(*weight), exWeights.setPixelType("UINT8")); - } - } - - float maxWeight = 0; - // image size - const int origWidth = weights[0]->width(); - const int origHeight = weights[0]->height(); - - // if we doing scaling, we have to backup L*a*b images of original size - std::vector<FLabImagePtr> backupLab; - if (flags & ADV_MULTIRES) { - for (unsigned int i = 0; i < LabImages.size(); i++) { - // backup original size L*a*b - backupLab.push_back(LabImages[i]); - } - } - - cout << endl << "Running khan algorithm" << endl; - // and we can run khan algorithm - // khan iteration - for (int it = 0; it < iterations; it++) { - if (verbosity > 0) - cout << "iteration " << it+1 << endl; - // copy weights from previous iteration - if (verbosity > 1) - cout << "copying weights from previous iteration" << endl; - - std::vector<FImagePtr> prevWeights; - for (unsigned int i = 0; i < weights.size(); i++) { - // scale weights to the requied size - if (flags & ADV_MULTIRES) { - // it would be better to use resampleImage, but it seems to be present only in VIGRA 1.6 - // so let's use some of the resizeImageINTERPOLATION() functions - - // compute width - int resized_width = origWidth / ( iterations/(it+1) ); - //compute height - int resized_height = origHeight / ( iterations/(it+1) ); - // destination images - FImage resizedWeight; - FLabImage resizedLab; - // it's not worthy to scale to less than 100px per side - if (resized_width > 100 && resized_height > 100) { - // create destination image of desired size - resizedWeight = FImage(Size2D(resized_width,resized_height)); - resizedLab = FLabImage(Size2D(resized_width,resized_height)); - } else if (origWidth >= 100 && origHeight >= 100) { - // resize it to the smallest value (ie 100px for the shorter side) - if (origWidth >= origHeight) { - resizedWeight = FImage(Size2D(100*origWidth/origHeight, 100)); - resizedLab = FLabImage(Size2D(100*origWidth/origHeight, 100)); - } else { - resizedWeight = FImage(Size2D(100, 100*origHeight/origWidth)); - resizedLab = FLabImage(Size2D(100, 100*origHeight/origWidth)); - } - } else { - // don't scale at all - // just copy weights as if no scaling seting was applied - goto DONTSCALE; - } - - // No interpolation – only for testing - resizeImageNoInterpolation(srcImageRange(*weights[i]), destImageRange(resizedWeight)); - resizeImageNoInterpolation(srcImageRange(*backupLab[i]), destImageRange(resizedLab)); - - FImagePtr tmp(new FImage(resizedWeight)); - prevWeights.push_back(tmp); - LabImages[i] = FLabImagePtr(new FLabImage(resizedLab)); - weights[i] = FImagePtr(new FImage(resizedWeight)); - } else { - DONTSCALE: - FImagePtr tmp(new FImage(*weights[i])); - prevWeights.push_back(tmp); - } - } - - // loop through all images - for (unsigned int i = 0; i < LabImages.size(); i++) { - if (verbosity > 1) - cout << "processing image " << i+1 << endl; - - // vector storing pixel data - AlgTinyVector<float, 3> X; - // sums for eq. 6 - double wpqssum = 0; - double wpqsKhsum = 0; - // image size - const int width = LabImages[i]->width(); - const int height = LabImages[i]->height(); - - // iterator to the upper left corner - FLabImage::traverser sy = LabImages[i]->upperLeft(); - // iterator to the lower right corner - FLabImage::traverser send = LabImages[i]->lowerRight(); - // iterator to the weight image left corner - FImage::traverser wy = weights[i]->upperLeft(); - // loop through the row - for (int y=0; sy.y != send.y; ++sy.y, ++wy.y, ++y) { - // iterator to the source (L*a*b image) - FLabImage::traverser sx = sy; - // iterator to the weight - FImage::traverser wx = wy; - // loop over the pixels - for (int x=0; sx.x != send.x; ++sx.x, ++wx.x, ++x) { - if (verbosity > 2) - cout << "processing pixel (" << x+1 << "," << y+1 << ")" << endl; - // set pixel vector - X = *sx; - - // loop through all layers - for (unsigned int j = 0; j < LabImages.size(); j++) { - if (verbosity > 2) - cout << "processing layer " << j << endl; - - // iterator to the neighbour - FLabImage::traverser neighby = LabImages[j]->upperLeft(); - // iterator to the weight - FImage::traverser weighty = prevWeights[j]->upperLeft(); - // pixel offset - int ndy = -NEIGHB_DIST; - // move them to the upper bound - // find the best upper bound - if (y-NEIGHB_DIST < 0) { - ndy = -y; - } - else { - neighby.y += y-NEIGHB_DIST; - weighty.y += y-NEIGHB_DIST; - } - - // iterate through neighbourhoods y axis - int maxDisty = (height - y) > NEIGHB_DIST ? NEIGHB_DIST : (height - y-1); - for (; ndy <= maxDisty; ++neighby.y, ++weighty.y, ++ndy) { - FLabImage::traverser neighbx = neighby; - FImage::traverser weightx = weighty; - // pixel offset - int ndx = -NEIGHB_DIST; - // move them to the upper bound - // find the best upper bound - if (x-NEIGHB_DIST < 0) { - ndx = -x; - } - else { - neighbx.x += x-NEIGHB_DIST; - weightx.x += x-NEIGHB_DIST; - } - // iterate through neighbourhoods x axis - int maxDistx = (width - x) > NEIGHB_DIST ? NEIGHB_DIST : (width - x-1); - for (; ndx <= maxDistx; ++neighbx.x, ++weightx.x, ++ndx) { - if (verbosity > 3) - cout << "(" << ndx << "," << ndy << ")"; - // now we can construct pixel vector - // should omit the middle pixel, ie use only neighbours - if (ndx != 0 || ndy != 0) { - wpqsKhsum += (*weightx * Kh(X-(*neighbx))); - wpqssum += *weightx; - } - - maxDistx = (width - x) > NEIGHB_DIST ? NEIGHB_DIST : (width - x-1); - } - if (verbosity > 3) - cout << endl; - - maxDisty = (height - y) > NEIGHB_DIST ? NEIGHB_DIST : (height - y-1); - } - } - - if (verbosity > 2) - cout << "computing new weight" << endl; - // compute probability and set weight - //cout << "P=" << (float) wpqsKhsum/wpqssum << endl; - if (flags & ADV_ONLYP) - *wx = (float) wpqsKhsum/wpqssum; - else - *wx *= (float) wpqsKhsum/wpqssum; - if (maxWeight < *wx) - maxWeight = *wx; - wpqsKhsum = wpqssum = 0; - - } - } - } - } - - if (verbosity > 1) - cout << "normalizing weights" << endl; - double factor = 255.0f/maxWeight; - for (unsigned int i=0; i<weights.size(); ++i) { - transformImage(srcImageRange(*(weights[i])), destImage(*(weights[i])), NormalizeFunctor<float>(factor)); - } - return weights; - } -} Modified: hugin/branches/gsoc2009_deghosting/src/deghosting/khan.h =================================================================== --- hugin/branches/gsoc2009_deghosting/src/deghosting/khan.h 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/deghosting/khan.h 2009-09-13 17:02:46 UTC (rev 4403) @@ -31,6 +31,9 @@ #include <vigra/transformimage.hxx> #include <vigra/colorconversions.hxx> +// for resampleImage +#include <vigra/resizeimage.hxx> + // for RGBvalue used in hat function #include <vigra/rgbvalue.hxx> @@ -54,14 +57,34 @@ namespace deghosting { - #ifdef DEGHOSTING_CACHE_IMAGES - typedef CachedFileImage<AlgTinyVector<float, 3> > FLabImage; - #else - typedef BasicImage<AlgTinyVector<float, 3> > FLabImage; - #endif - typedef boost::shared_ptr<FLabImage> FLabImagePtr; - - class Khan : public Deghosting + template <class PixelType> + class ImageTypes { + public: + typedef vigra::FImage ImageType; + typedef boost::shared_ptr<ImageType> ImagePtr; + #ifdef DEGHOSTING_CACHE_IMAGES + typedef CachedFileImage<float> ProcessImageType; + #else + typedef BasicImage<float> ProcessImageType; + #endif + typedef boost::shared_ptr<ProcessImageType> ProcessImageTypePtr; + }; + + template <class PixelType> + class ImageTypes<RGBValue<PixelType> > { + public: + typedef vigra::FRGBImage ImageType; + typedef boost::shared_ptr<ImageType> ImagePtr; + #ifdef DEGHOSTING_CACHE_IMAGES + typedef CachedFileImage<AlgTinyVector<float, 3> > ProcessImageType; + #else + typedef BasicImage<AlgTinyVector<float, 3> > ProcessImageType; + #endif + typedef boost::shared_ptr<ProcessImageType> ProcessImageTypePtr; + }; + + template <class PixelType> + class Khan : public Deghosting, private ImageTypes<PixelType> { public: Khan(std::vector<std::string>& inputFiles, const uint16_t flags, const uint16_t debugFlags, int iterations, double sigma, int verbosity); @@ -69,6 +92,12 @@ std::vector<FImagePtr> createWeightMasks(); ~Khan() {} protected: + typedef typename ImageTypes<PixelType>::ProcessImageType ProcessImageType; + typedef typename ImageTypes<PixelType>::ProcessImageType::traverser ProcessImageTraverser; + typedef typename ImageTypes<PixelType>::ProcessImageType::PixelType ProcessImagePixelType; + typedef typename ImageTypes<PixelType>::ProcessImageTypePtr ProcessImageTypePtr; + typedef typename NumericTraits<PixelType>::isScalar srcIsScalar; + // Kh() things // (2*pi)^(1/2) double PIPOW; @@ -78,7 +107,7 @@ double sigma; // other necessary stuff - std::vector<FLabImagePtr> LabImages; + std::vector<ProcessImageTypePtr> LabImages; std::vector<FImagePtr> weights; /** set sigma @@ -95,15 +124,377 @@ /** kernel function * Standard probability density function */ - inline float Kh(deghosting::AlgTinyVector< float, 3 > x); + inline float Kh(ProcessImagePixelType x); /** function to preprocess input image * This function loads image, linearize it using EMoR (FIXME), * tranform it using logarithm or gamma if input images are HDR */ - template <class PixelType> - void preprocessImage(unsigned int i, FImagePtr &weight, FLabImagePtr &LabImage); + void preprocessImage(unsigned int i, FImagePtr &weight, ProcessImageTypePtr &LabImage); }; + + template <class PixelType> + Khan<PixelType>::Khan(std::vector<std::string>& newInputFiles, const uint16_t newFlags, const uint16_t newDebugFlags, + int newIterations, double newSigma, int newVerbosity) { + try { + Deghosting::loadImages(newInputFiles); + Deghosting::setFlags(newFlags); + Deghosting::setDebugFlags(newDebugFlags); + Deghosting::setIterationNum(newIterations); + Deghosting::setVerbosity(newVerbosity); + + // I don't know why, but sigma for HDR input have to approximately 10 times smaller + // FIXME: Maybe it would be better to use different sigma for different images in case both HDR and LDR are mixed + const char * fileType= ImageImportInfo(newInputFiles[0].c_str()).getFileType(); + if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { + setSigma(newSigma/10); + } else { + setSigma(newSigma); + } + + for (unsigned int i=0; i<5; i++) + Deghosting::response.push_back(0); + PIPOW = sigma*std::sqrt(2*PI); + denom = 1/PIPOW; + } catch (...) { + throw; + } + } + + template <class PixelType> + Khan<PixelType>::Khan(std::vector<ImageImportInfo>& newInputFiles, const uint16_t newFlags, const uint16_t newDebugFlags, + int newIterations, double newSigma, int newVerbosity) { + try { + Deghosting::loadImages(newInputFiles); + Deghosting::setFlags(newFlags); + Deghosting::setDebugFlags(newDebugFlags); + Deghosting::setIterationNum(newIterations); + Deghosting::setVerbosity(newVerbosity); + + // I don't know why, but sigma for HDR input have to approximately 10 times smaller + // FIXME: Maybe it would be better to use different sigma for different images in case both HDR and LDR are mixed + const char * fileType= newInputFiles[0].getFileType(); + if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { + setSigma(newSigma/10); + } else { + setSigma(newSigma); + } + + for (unsigned int i=0; i<5; i++) + Deghosting::response.push_back(0); + PIPOW = sigma*std::sqrt(2*PI); + denom = 1/PIPOW; + } catch (...) { + throw; + } + } + + template <class PixelType> + void Khan<PixelType>::setSigma(double newSigma) { + sigma = newSigma; + } + + template <class PixelType> + float Khan<PixelType>::Kh(ProcessImagePixelType x) { + #ifdef ATAN_KH + // good choice for sigma for this function is around 600 + return std::atan(-(x*x)+sigma)/PI + 0.5; + #else + // good choice for sigma for this function is around 30 + return (std::exp(-(x*x)/(2*sigma*sigma)) * denom); + #endif + } + + /*void Khan::linearizeRGB(std::string inputFile,FRGBImage *pInputImg) { + HuginBase::SrcPanoImage panoImg(inputFile); + panoImg.setResponseType(HuginBase::SrcPanoImage::RESPONSE_EMOR); + panoImg.setEMoRParams(response); + // response transform functor + HuginBase::Photometric::InvResponseTransform<RGBValue<float>, + RGBValue<float> > invResponse(panoImg); + invResponse.enforceMonotonicity(); + + // iterator to the upper left corner + FRGBImage::traverser imgIterSourceY = pInputImg->upperLeft(); + // iterator to he lower right corner + FRGBImage::traverser imgIterEnd = pInputImg->lowerRight(); + // iterator to the output + FRGBImage::traverser imgIterOut = pInputImg->upperLeft(); + // loop through the image + for (int y=1; imgIterSourceY.y != imgIterEnd.y; ++imgIterSourceY.y, ++imgIterOut.y, ++y) { + // iterator to the input + FRGBImage::traverser sx = imgIterSourceY; + // iterator to the output + FRGBImage::traverser dx = imgIterOut; + for (int x=1; sx.x != imgIterEnd.x; ++sx.x, ++dx.x, ++x) { + // transform image using response + *dx = vigra_ext::zeroNegative(invResponse(*sx, hugin_utils::FDiff2D(x, y))); + } + } + }*/ + + template <class PixelType> + void Khan<PixelType>::preprocessImage(unsigned int i, FImagePtr &weight, ProcessImageTypePtr &LabImage) { + ImageImportInfo imgInfo(inputFiles[i]); + BasicImage<PixelType> * pInputImg = new BasicImage<PixelType>(imgInfo.size()); + BImage imgAlpha(imgInfo.size()); + weight = FImagePtr(new FImage(imgInfo.size())); + + // import alpha + importImageAlpha(imgInfo, destImage(*pInputImg),destImage(imgAlpha)); + + // linearize RGB and convert it to L*a*b image + //linearizeRGB(inputFiles[i], pInputImg); + + // take logarithm or gamma correction if the input images are HDR + // I'm not sure if it's the right way how to + // find out if they are HDR + const char * fileType= imgInfo.getFileType(); + if ( (!strcmp(fileType,"TIFF") && strcmp(fileType,"UINT8")) || !strcmp(fileType,"EXR") || !strcmp(fileType,"FLOAT")) { + // use gamma 2.2 + if (flags & ADV_GAMMA) { + // GammaFunctor is only in vigra 1.6 GRRR + // I have to use BrightnessContrastFunctor + // TODO: change to the GammaFunctor in the future + vigra::FindMinMax<float> minmax; + vigra::inspectImage(srcImageRange(*pInputImg), minmax); + transformImage(srcImageRange(*pInputImg),destImage(*pInputImg),BrightnessContrastFunctor<PixelType>(0.45,1.0,minmax.min, minmax.max)); + } else { + // take logarithm + transformImage(srcImageRange(*pInputImg),destImage(*pInputImg),LogarithmFunctor<PixelType>(1.0)); + } + } + + // generate initial weights using hat function + if (flags & ADV_NOINITWEIGHTS) { + // convert image to grayscale (and take logarithm) + RGBToGrayAccessor<FRGBImage::PixelType> color2gray; + transformImage(srcImageRange(*pInputImg, color2gray), destImage(*weight), log(Arg1()+Param(1.0f))); + } + else { + transformImage(srcImageRange(*pInputImg),destImage(*weight),HatFunctor<PixelType>()); + } + + //if (srcIsScalar()) { + // copyImage(srcImageRange(*pInputImg), destImage(*LabImage)); + //} else { + // convert from linear RGB to L*a*b // + RGB2LabFunctor<float> RGB2Lab; + LabImage = ProcessImageTypePtr(new ProcessImageType(imgInfo.size())); + transformImage(srcImageRange(*pInputImg), destImage(*LabImage), RGB2Lab); + //} + + delete pInputImg; + pInputImg = 0; + } + + template <class PixelType> + std::vector<FImagePtr> Khan<PixelType>::createWeightMasks() { + for (unsigned int i = 0; i < inputFiles.size(); i++) { + FImagePtr weight; + ProcessImageTypePtr LabImage; + preprocessImage(i, weight, LabImage); + LabImages.push_back(LabImage); + weights.push_back(weight); + + // save init weights + if (debugFlags & SAVE_INITWEIGHTS) { + char tmpfn[100]; + snprintf(tmpfn, 99, "init_weights_%d.tiff", i); + ImageExportInfo exWeights(tmpfn); + exportImage(srcImageRange(*weight), exWeights.setPixelType("UINT8")); + } + } + + float maxWeight = 0; + // image size + const int origWidth = weights[0]->width(); + const int origHeight = weights[0]->height(); + + // if we doing scaling, we have to backup L*a*b images of original size + std::vector<ProcessImageTypePtr> backupLab; + if (flags & ADV_MULTIRES) { + for (unsigned int i = 0; i < LabImages.size(); i++) { + // backup original size L*a*b + backupLab.push_back(LabImages[i]); + } + } + + cout << endl << "Running khan algorithm" << endl; + // and we can run khan algorithm + // khan iteration + for (int it = 0; it < iterations; it++) { + if (verbosity > 0) + cout << "iteration " << it+1 << endl; + // copy weights from previous iteration + if (verbosity > 1) + cout << "copying weights from previous iteration" << endl; + + std::vector<FImagePtr> prevWeights; + for (unsigned int i = 0; i < weights.size(); i++) { + // scale weights to the requied size + if (flags & ADV_MULTIRES) { + // it would be better to use resampleImage, but it seems to be present only in VIGRA 1.6 + // so let's use some of the resizeImageINTERPOLATION() functions + + // compute width + int resized_width = origWidth / ( iterations/(it+1) ); + //compute height + int resized_height = origHeight / ( iterations/(it+1) ); + // destination images + FImage resizedWeight; + ProcessImageType resizedLab; + // it's not worthy to scale to less than 100px per side + if (resized_width > 100 && resized_height > 100) { + // create destination image of desired size + resizedWeight = FImage(Size2D(resized_width,resized_height)); + resizedLab = ProcessImageType(Size2D(resized_width,resized_height)); + } else if (origWidth >= 100 && origHeight >= 100) { + // resize it to the smallest value (ie 100px for the shorter side) + if (origWidth >= origHeight) { + resizedWeight = FImage(Size2D(100*origWidth/origHeight, 100)); + resizedLab = ProcessImageType(Size2D(100*origWidth/origHeight, 100)); + } else { + resizedWeight = FImage(Size2D(100, 100*origHeight/origWidth)); + resizedLab = ProcessImageType(Size2D(100, 100*origHeight/origWidth)); + } + } else { + // don't scale at all + // just copy weights as if no scaling seting was applied + goto DONTSCALE; + } + + // No interpolation – only for testing + resizeImageNoInterpolation(srcImageRange(*weights[i]), destImageRange(resizedWeight)); + resizeImageNoInterpolation(srcImageRange(*backupLab[i]), destImageRange(resizedLab)); + + FImagePtr tmp(new FImage(resizedWeight)); + prevWeights.push_back(tmp); + LabImages[i] = ProcessImageTypePtr(new ProcessImageType(resizedLab)); + weights[i] = FImagePtr(new FImage(resizedWeight)); + } else { + DONTSCALE: + FImagePtr tmp(new FImage(*weights[i])); + prevWeights.push_back(tmp); + } + } + + // loop through all images + for (unsigned int i = 0; i < LabImages.size(); i++) { + if (verbosity > 1) + cout << "processing image " << i+1 << endl; + + // vector storing pixel data + AlgTinyVector<float, 3> X; + // sums for eq. 6 + double wpqssum = 0; + double wpqsKhsum = 0; + // image size + const int width = LabImages[i]->width(); + const int height = LabImages[i]->height(); + + // iterator to the upper left corner + ProcessImageTraverser sy = LabImages[i]->upperLeft(); + // iterator to the lower right corner + ProcessImageTraverser send = LabImages[i]->lowerRight(); + // iterator to the weight image left corner + FImage::traverser wy = weights[i]->upperLeft(); + // loop through the row + for (int y=0; sy.y != send.y; ++sy.y, ++wy.y, ++y) { + // iterator to the source (L*a*b image) + ProcessImageTraverser sx = sy; + // iterator to the weight + FImage::traverser wx = wy; + // loop over the pixels + for (int x=0; sx.x != send.x; ++sx.x, ++wx.x, ++x) { + if (verbosity > 2) + cout << "processing pixel (" << x+1 << "," << y+1 << ")" << endl; + // set pixel vector + X = *sx; + + // loop through all layers + for (unsigned int j = 0; j < LabImages.size(); j++) { + if (verbosity > 2) + cout << "processing layer " << j << endl; + + // iterator to the neighbour + ProcessImageTraverser neighby = LabImages[j]->upperLeft(); + // iterator to the weight + FImage::traverser weighty = prevWeights[j]->upperLeft(); + // pixel offset + int ndy = -NEIGHB_DIST; + // move them to the upper bound + // find the best upper bound + if (y-NEIGHB_DIST < 0) { + ndy = -y; + } + else { + neighby.y += y-NEIGHB_DIST; + weighty.y += y-NEIGHB_DIST; + } + + // iterate through neighbourhoods y axis + int maxDisty = (height - y) > NEIGHB_DIST ? NEIGHB_DIST : (height - y-1); + for (; ndy <= maxDisty; ++neighby.y, ++weighty.y, ++ndy) { + ProcessImageTraverser neighbx = neighby; + FImage::traverser weightx = weighty; + // pixel offset + int ndx = -NEIGHB_DIST; + // move them to the upper bound + // find the best upper bound + if (x-NEIGHB_DIST < 0) { + ndx = -x; + } + else { + neighbx.x += x-NEIGHB_DIST; + weightx.x += x-NEIGHB_DIST; + } + // iterate through neighbourhoods x axis + int maxDistx = (width - x) > NEIGHB_DIST ? NEIGHB_DIST : (width - x-1); + for (; ndx <= maxDistx; ++neighbx.x, ++weightx.x, ++ndx) { + if (verbosity > 3) + cout << "(" << ndx << "," << ndy << ")"; + // now we can construct pixel vector + // should omit the middle pixel, ie use only neighbours + if (ndx != 0 || ndy != 0) { + wpqsKhsum += (*weightx * Kh(X-(*neighbx))); + wpqssum += *weightx; + } + + maxDistx = (width - x) > NEIGHB_DIST ? NEIGHB_DIST : (width - x-1); + } + if (verbosity > 3) + cout << endl; + + maxDisty = (height - y) > NEIGHB_DIST ? NEIGHB_DIST : (height - y-1); + } + } + + if (verbosity > 2) + cout << "computing new weight" << endl; + // compute probability and set weight + //cout << "P=" << (float) wpqsKhsum/wpqssum << endl; + if (flags & ADV_ONLYP) + *wx = (float) wpqsKhsum/wpqssum; + else + *wx *= (float) wpqsKhsum/wpqssum; + if (maxWeight < *wx) + maxWeight = *wx; + wpqsKhsum = wpqssum = 0; + + } + } + } + } + + if (verbosity > 1) + cout << "normalizing weights" << endl; + double factor = 255.0f/maxWeight; + for (unsigned int i=0; i<weights.size(); ++i) { + transformImage(srcImageRange(*(weights[i])), destImage(*(weights[i])), NormalizeFunctor<float>(factor)); + } + return weights; + } } #endif /* KHAN_H_ */ Modified: hugin/branches/gsoc2009_deghosting/src/tools/CMakeLists.txt =================================================================== --- hugin/branches/gsoc2009_deghosting/src/tools/CMakeLists.txt 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/tools/CMakeLists.txt 2009-09-13 17:02:46 UTC (rev 4403) @@ -23,7 +23,7 @@ add_executable(pto2mk pto2mk.cpp) target_link_libraries(pto2mk ${common_libs} ${image_libs} ) -add_executable(hugin_hdrmerge hugin_hdrmerge.cpp ../deghosting/deghosting.cpp ../deghosting/khan.cpp) +add_executable(hugin_hdrmerge hugin_hdrmerge.cpp ../deghosting/deghosting.cpp) target_link_libraries(hugin_hdrmerge ${common_libs} ${image_libs} ) install(TARGETS nona vig_optimize autooptimiser fulla align_image_stack tca_correct pto2mk hugin_hdrmerge Modified: hugin/branches/gsoc2009_deghosting/src/tools/hugin_hdrmerge.cpp =================================================================== --- hugin/branches/gsoc2009_deghosting/src/tools/hugin_hdrmerge.cpp 2009-09-13 15:33:18 UTC (rev 4402) +++ hugin/branches/gsoc2009_deghosting/src/tools/hugin_hdrmerge.cpp 2009-09-13 17:02:46 UTC (rev 4403) @@ -293,7 +293,7 @@ reduceFilesToHDR(inputFiles, outputFile, onlyCompleteOverlap, waverage); } else if (mode == "khan") { Deghosting* deghoster = NULL; - Khan khanDeghoster(inputFiles, 0, 0, iterations, sigma, g_verbose); + Khan<RGBValue<float> > khanDeghoster(inputFiles, 0, 0, iterations, sigma, g_verbose); deghoster = &khanDeghoster; vector<FImagePtr> weights = deghoster->createWeightMasks(); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |