NearestFeatureTransform.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00002 /*
00003  * Copyright (C) 2004 Andrew Mihal
00004  *
00005  * This file is part of Enblend.
00006  *
00007  * Adapted to vigra & panoramic stitcher by Pablo d'Angelo
00008  *
00009  * Enblend is free software; you can redistribute it and/or modify
00010  * it under the terms of the GNU General Public License as published by
00011  * the Free Software Foundation; either version 2 of the License, or
00012  * (at your option) any later version.
00013  *
00014  * Enblend is distributed in the hope that it will be useful,
00015  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017  * GNU General Public License for more details.
00018  *
00019  * You should have received a copy of the GNU General Public License
00020  * along with Enblend; if not, write to the Free Software
00021  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022  */
00023 
00024 #ifndef NEARESTFEATURETRANSFORM_H
00025 #define NEARESTFEATURETRANSFORM_H
00026 
00027 #include <iostream>
00028 #include <list>
00029 #include <math.h>
00030 #include <vigra/basicimage.hxx>
00031 #include <vigra/functorexpression.hxx>
00032 #include <vigra/combineimages.hxx>
00033 
00034 #include <hugin_utils/utils.h>
00035 #include <appbase/ProgressDisplayOld.h>
00036 
00037 // HACK, for uint typedefs...
00038 #include <tiff.h>
00039 #include <tiffio.h>
00040 
00041 //#include "enblend.h"
00042 
00043 
00044 
00045 /*
00046 extern int Verbose;
00047 extern bool Wraparound;
00048 extern uint32 OutputWidth;
00049 extern uint32 OutputHeight;
00050 
00051 // Region of interest for this operation.
00052 extern uint32 UBBFirstX;
00053 extern uint32 UBBLastX;
00054 extern uint32 UBBFirstY;
00055 extern uint32 UBBLastY;
00056 */
00057 
00058 namespace vigra_ext {
00059 
00060 #define EUCLIDEAN_METRIC
00061 
00062 #ifdef EUCLIDEAN_METRIC
00063 //    typedef uint32 dist_t;
00064 //    #define DIST_MAX ((dist_t)-1)
00065     typedef float dist_t;
00066     #define DIST_MAX (1e20)
00067     #define DIST_MIN 0
00068     inline dist_t distance(uint32 deltaX, dist_t *deltaY) {
00069         return (*deltaY == DIST_MAX) ? DIST_MAX : *deltaY + deltaX * deltaX;
00070     }
00071     inline dist_t distance(uint32 deltaY) {
00072         return deltaY * deltaY;
00073     }
00074 #else
00075 #ifdef CHESSBOARD_METRIC
00076     typedef uint32 dist_t;
00077     #define DIST_MAX ((dist_t)-1)
00078     #define DIST_MIN 0
00079     inline dist_t distance(uint32 deltaX, dist_t *deltaY) {
00080         return max(deltaX, *deltaY);
00081     }
00082     inline dist_t distance(uint32 deltaY) {
00083         return deltaY;
00084     }
00085 #else
00086 #ifdef MANHATTAN_METRIC
00087     typedef uint32 dist_t;
00088     #define DIST_MAX ((dist_t)-1)
00089     #define DIST_MIN 0
00090     inline dist_t distance(uint32 deltaX, dist_t *deltaY) {
00091         return (*deltaY == DIST_MAX) ? DIST_MAX : deltaX + *deltaY;
00092     }
00093     inline dist_t distance(uint32 deltaY) {
00094         return deltaY;
00095     }
00096 #endif
00097 #endif
00098 #endif
00099 
00100 
00101 #ifdef DEBUG_NEAREST
00102 static inline void saveDist(dist_t *ptr, uint32 w, uint32 h, const char * filename)
00103 {
00104 
00105     TIFF *outputTIFF = TIFFOpen(filename, "w");
00106     if (outputTIFF == NULL) {
00107         std::cerr << "enblend: error opening TIFF file \""
00108              << filename << "\"" << std::endl;
00109         exit(1);
00110     }
00111 
00112     TIFFSetField(outputTIFF, TIFFTAG_ORIENTATION, 1);
00113     TIFFSetField(outputTIFF, TIFFTAG_SAMPLEFORMAT, 3);
00114     TIFFSetField(outputTIFF, TIFFTAG_SAMPLESPERPIXEL, 1);
00115     TIFFSetField(outputTIFF, TIFFTAG_BITSPERSAMPLE, 32);
00116     TIFFSetField(outputTIFF, TIFFTAG_IMAGEWIDTH, w);
00117     TIFFSetField(outputTIFF, TIFFTAG_IMAGELENGTH, h);
00118     TIFFSetField(outputTIFF, TIFFTAG_PHOTOMETRIC, 1);
00119     TIFFSetField(outputTIFF, TIFFTAG_PLANARCONFIG, 1);
00120 
00121     for (uint32 scan = 0; scan < h; scan++) {
00122         TIFFWriteScanline(outputTIFF,
00123                           &(ptr[scan * w]),
00124                           scan,
00125                           1);
00126     }
00127     TIFFClose(outputTIFF);
00128 }
00129 
00130 static void saveDistRaw(void * ptr, uint32 w, uint32 h, const char * file)
00131 {
00132     FILE * f = fopen(file,"wb+");
00133     if (f == NULL) {
00134         std::cerr << "enblend: error creating temporary file." << std::endl;
00135         exit(1);
00136     }
00137     size_t size = w*h;
00138     size_t itemsWritten = fwrite(ptr, sizeof(dist_t), size, f);
00139     if (itemsWritten < size) {
00140         std::cerr << "enblend: error writing to temporary file." << std::endl;
00141         perror("reason");
00142         exit(1);
00143     }
00144 }
00145 #endif
00146 
00155 template <typename Feat1Iter, typename Feat1Accessor,
00156           typename Feat2Iter, typename Feat2Accessor,
00157           typename Mask1Iter, typename Mask1Accessor,
00158           typename Mask2Iter, typename Mask2Accessor>
00159 void nearestFeatureTransform(vigra::triple<Feat1Iter, Feat1Iter, Feat1Accessor> feat1,
00160                              std::pair<Feat2Iter, Feat2Accessor> feat2,
00161                              std::pair<Mask1Iter, Mask1Accessor> mask1,
00162                              std::pair<Mask2Iter, Mask2Accessor> mask2,
00163                              AppBase::MultiProgressDisplay & progress)
00164 {
00165     typedef typename Mask1Accessor::value_type MaskType;
00166     typedef vigra::BasicImage<MaskType>       FeatImage;
00167     typedef typename FeatImage::Iterator      FeatIter;
00168     progress.pushTask(AppBase::ProgressTask("blend mask", "creating blend line",1/4));  
00169 
00170     vigra::Diff2D ubbSize = feat1.second - feat1.first;
00171 
00172     uint32 ubbWidth = ubbSize.x;
00173     uint32 ubbHeight = ubbSize.y;
00174     uint32 ubbPixels = ubbWidth * ubbHeight;
00175 
00176     // mask for the feature pixels ( a feature pixel is
00177     // feat1 xor feat2
00178     FeatImage feature(ubbSize);
00179 
00180     // For each pixel, store the distance to the nearest feature in the same
00181     // column.
00182     dist_t *dnfColumn = (dist_t*)malloc(ubbPixels * sizeof(dist_t));
00183 
00184     if (dnfColumn == NULL) {
00185         std::cerr << std::endl
00186              << "enblend: out of memory (in nearestFeatureTransform for dnfColumn)" << std::endl;
00187         exit(1);
00188     }
00189 
00190     // For each pixel, store the distance to the nearest feature in the same
00191     // column or any column to the left.
00192     dist_t *dnfLeft = (dist_t*)malloc(ubbPixels * sizeof(dist_t));
00193     if (dnfLeft == NULL) {
00194         std::cerr << "enblend: out of memory (in nearestFeatureTransform for dnfLeft)" << std::endl;
00195         exit(1);
00196     }
00197 
00198     // Initialize dnfColumn top-down. Store the distance to the nearest feature
00199     // in the same column and above us.
00200     /*
00201     if (Verbose > 0) {
00202         cout << "Creating blend mask: 1/4 ";
00203         cout.flush();
00204     }
00205     */
00206 
00207     //    MaskPixel *firstMaskP = mask;
00208 
00209     Mask1Iter m1x(mask1.first);
00210     Mask2Iter m2x(mask2.first);
00211 
00212     Feat1Iter f1x(feat1.first);
00213     Feat2Iter f2x(feat2.first);
00214 
00215     //FeatIter featx(feature.upperLeft());
00216     //    MaskIter mux(maskUnion);
00217 
00218     for (uint32 x = 0; x < ubbWidth; x++, ++m1x.x, ++m2x.x, ++f1x.x, ++f2x.x) {
00219         dist_t *dnfColumnP = &dnfColumn[x];
00220         //MaskPixel *maskP = firstMaskP + x;
00221 
00222         // Color of the last feature pixel.
00223         MaskType lastFeature2 = 0;
00224         MaskType lastFeature1 = 0;
00225         // Distance to the last feature pixel in pixels.
00226         uint32 lastFeatureDeltaY = 0;
00227         bool foundFirstFeature = false;
00228 
00229         Mask1Iter m1y(m1x);
00230         Mask2Iter m2y(m2x);
00231         Feat1Iter f1y(f1x);
00232         Feat2Iter f2y(f2x);
00233         for (uint32 y = 0; y < ubbHeight; y++, ++m1y.y, ++m2y.y, ++f1y.y, ++f2y.y) {
00234             // initialize masks and feature mask (pixel defined by only one image).
00235             if ((*f2y) != (*f1y)) {
00236             //if (maskP->r == 0) {
00237                 // maskP is a feature pixel
00238                 *dnfColumnP = DIST_MIN;
00239                 lastFeatureDeltaY = 0;
00240                 //lastFeatureG = maskP->g;
00241                 //lastFeatureB = maskP->b;
00242                 lastFeature2 = *f2y;
00243                 lastFeature1 = *f1y;
00244                 *m1y = lastFeature1;
00245                 *m2y = lastFeature2;
00246                 foundFirstFeature = true;
00247             } else if (foundFirstFeature) {
00248                 // maskP is not a feature.
00249                 *dnfColumnP = distance(lastFeatureDeltaY);
00250                 //maskP->g = lastFeatureG;
00251                 //maskP->b = lastFeatureB;
00252                 *m1y = lastFeature1;
00253                 *m2y = lastFeature2;
00254             } else {
00255                 *dnfColumnP = DIST_MAX;
00256             }
00257 
00258             lastFeatureDeltaY++;
00259 
00260             // Move pointers down one row.
00261             dnfColumnP += ubbWidth;
00262             //            maskP += ubbWidth;
00263         }
00264     }
00265 
00266 #ifdef DEBUG_NEAREST
00267     vigra::exportImage(vigra::srcImageRange(feature), vigra::ImageExportInfo("nona_0_featuremask.tif"));
00268     vigra::exportImage(srcIterRange(mask1.first, mask1.first + ubbSize), vigra::ImageExportInfo("nona_1_mask1.tif"));
00269     vigra::exportImage(srcIterRange(mask2.first, mask2.first + ubbSize), vigra::ImageExportInfo("nona_1_mask2.tif"));
00270     saveDist(dnfColumn, ubbWidth, ubbHeight, "nona_1_dnfColumn.tif");
00271     saveDist(dnfLeft, ubbWidth, ubbHeight, "nona_1_dnfLeft.tif");
00272 #endif
00273 
00274     // Initialize dnfColumn bottom-up. Caluclate the distance to the nearest
00275     // feature in the same column and below us.
00276     // If this is smaller than the value caluclated in the top-down pass,
00277     // overwrite that value.
00278     progress.increase();
00279     /*
00280     if (Verbose > 0) {
00281         cout << "2/4 ";
00282         cout.flush();
00283     }
00284     */
00285 
00286     //MaskPixel *lastMaskP = &mask[ubbPixels - 1];
00287 
00288     m1x = mask1.first + ubbSize - vigra::Diff2D(1,1);
00289     m2x = mask2.first + ubbSize - vigra::Diff2D(1,1);
00290     f1x = feat1.first + ubbSize - vigra::Diff2D(1,1);
00291     f2x = feat2.first + ubbSize - vigra::Diff2D(1,1);
00292 
00293     dist_t *lastDNFColumnP = &dnfColumn[ubbPixels - 1];
00294     for (uint32 x = 0; x < ubbWidth; x++, --m1x.x, --m2x.x, --f1x.x, --f2x.x) {
00295         dist_t *dnfColumnP = lastDNFColumnP - x;
00296         //MaskPixel *maskP = lastMaskP - x;
00297 
00298         // Color of the last feature pixel.
00299         MaskType lastFeatureG = 0;
00300         MaskType lastFeatureB = 0;
00301         // Distance to the last feature pixel in pixels.
00302         uint32 lastFeatureDeltaY = 0;
00303         bool foundFirstFeature = false;
00304 
00305         Mask1Iter m1y(m1x);
00306         Mask2Iter m2y(m2x);
00307         Feat1Iter f1y(f1x);
00308         Feat2Iter f2y(f2x);
00309         for (uint32 y = 0; y < ubbHeight; y++, --m1y.y, --m2y.y, --f1y.y, --f2y.y)
00310         {
00311             //if (maskP->r == 0) {
00312             if ((*f2y) != (*f1y)) {
00313                 // maskP is a feature pixel.
00314                 //*dnfColumnP = DIST_MIN; don't need to do this again.
00315                 lastFeatureDeltaY = 0;
00316                 //lastFeatureG = maskP->g;
00317                 //lastFeatureB = maskP->b;
00318                 lastFeatureG = *m2y;
00319                 lastFeatureB = *m1y;
00320                 foundFirstFeature = true;
00321             } else if (foundFirstFeature) {
00322                 // maskP is not a feature.
00323                 dist_t distLastFeature = distance(lastFeatureDeltaY);
00324                 // If last feature is closer than nearest feature above,
00325                 // change distance and color to match last feature.
00326                 if (distLastFeature < *dnfColumnP) {
00327                     *dnfColumnP = distLastFeature;
00328                     //maskP->g = lastFeatureG;
00329                     //maskP->b = lastFeatureB;
00330                     *m1y = lastFeatureB;
00331                     *m2y = lastFeatureG;
00332                 }
00333             }
00334 
00335             lastFeatureDeltaY++;
00336 
00337             // Move pointers up one row.
00338             dnfColumnP -= ubbWidth;
00339             //maskP -= ubbWidth;
00340         }
00341     }
00342 
00343 #ifdef DEBUG_NEAREST
00344     vigra::exportImage(srcIterRange(mask1.first, mask1.first + ubbSize), vigra::ImageExportInfo("nona_2_mask1.tif"));
00345     vigra::exportImage(srcIterRange(mask2.first, mask2.first + ubbSize), vigra::ImageExportInfo("nona_2_mask2.tif"));
00346     saveDist(dnfColumn, ubbWidth, ubbHeight, "nona_2_dnfColumn.tif");
00347     saveDist(dnfLeft, ubbWidth, ubbHeight, "nona_2_dnfLeft.tif");
00348 #endif
00349     
00350     //size_t maxListSize = 0;
00351 
00352     // Calculate dnfLeft for each pixel.
00353     progress.increase();
00354     /*
00355     if (Verbose > 0) {
00356         cout << "3/4 ";
00357         cout.flush();
00358     }
00359     */
00360 
00361     Mask1Iter m1y(mask1.first);
00362     Mask2Iter m2y(mask2.first);
00363     for (uint32 y = 0; y < ubbHeight; y++, ++m1y.y, ++m2y.y) {
00364         dist_t *dnfLeftP = &dnfLeft[y * ubbWidth];
00365         dist_t *dnfColumnP = &dnfColumn[y * ubbWidth];
00366         //MaskPixel *maskP = firstMaskP + (y * ubbWidth);
00367 
00368         // List of dnfColumnP's on the left that might be the closest features
00369         // to the current dnfColumnP.
00370         std::list<dist_t*> potentialFeatureList;
00371 
00372         Mask1Iter m1x(m1y);
00373         Mask2Iter m2x(m2y);
00374         for (uint32 x = 0; x < ubbWidth; x++, ++m1x.x, ++m2x.x) {
00375             // First add ourself to the list.
00376             potentialFeatureList.push_back(dnfColumnP);
00377 
00378             // Iterate through the list starting at the right. For each
00379             // potential feature, all of the potential features to the left
00380             // in the list must be strictly closer. If not delete them from
00381             // the list.
00382             std::list<dist_t*>::iterator potentialFeature =
00383                     --(potentialFeatureList.end());
00384             // The last potential feature is dnfColumnP, just added above.
00385             // That is in the current column so the distance to that feature
00386             // is simply *dnfColumnP.
00387             dist_t distPotentialFeature = *dnfColumnP;
00388             while (potentialFeature != potentialFeatureList.begin()) {
00389                 // Make an iterator that points to the predecessor.
00390                 std::list<dist_t*>::iterator previousFeature = potentialFeature;
00391                 previousFeature--;
00392 
00393                 // previousFeature is this many columns to the left of (x,y).
00394                 uint32 deltaX = dnfColumnP - *previousFeature;
00395 
00396                 // previousFeature is this far from (x,y).
00397                 dist_t distPreviousFeature = distance(deltaX, *previousFeature);
00398 
00399                 if (distPreviousFeature >= distPotentialFeature) {
00400                     // previousFeature is not a candidate for dnfLeftP
00401                     // or anything further to the right.
00402                     potentialFeatureList.erase(previousFeature);
00403                 } else {
00404                     // previousFeature is a candidate.
00405                     potentialFeature = previousFeature;
00406                     distPotentialFeature = distPreviousFeature;
00407                 }
00408             }
00409 
00410             // The closest feature to (x,y) in columns <= x is the first
00411             // potential feature in the list.
00412             *dnfLeftP = distPotentialFeature;
00413 
00414             // Set color of maskP to be color of closest feature to the left.
00415             //MaskPixel *maskPLeft = maskP - (dnfColumnP - *potentialFeature);
00416             int dx = (dnfColumnP - *potentialFeature);
00417             Mask1Iter m1l(m1x);
00418             m1l.x -= dx;
00419             *m1x = *m1l;
00420             Mask1Iter m2l(m2x);
00421             m2l.x -= dx;
00422             *m2x = *m2l;
00423 
00424             // Move pointers right one column.
00425             dnfLeftP++;
00426             dnfColumnP++;
00427             //maskP++;
00428 
00429             //maxListSize = max(maxListSize, potentialFeatureList.size());
00430         }
00431     }
00432     //if (Verbose > 0) {
00433     //    cout << "max feature list size=" << maxListSize << endl;
00434     //}
00435     //maxListSize = 0;
00436 
00437 #ifdef DEBUG_NEAREST    
00438     vigra::exportImage(vigra::srcIterRange(mask1.first, mask1.first + ubbSize), vigra::ImageExportInfo("nona_3_mask1.tif"));
00439     vigra::exportImage(vigra::srcIterRange(mask2.first, mask2.first + ubbSize), vigra::ImageExportInfo("nona_3_mask2.tif"));
00440     saveDist(dnfColumn, ubbWidth, ubbHeight, "nona_3_dnfColumn.tif");
00441     saveDist(dnfLeft, ubbWidth, ubbHeight, "nona_3_dnfLeft.tif");
00442 #endif
00443     
00444     // Final pass: calculate the distance to the nearest feature in the same
00445     // column or any column to the right. If this is smaller than dnfLeftP,
00446     // Then recolor the pixel to the color of the nearest feature to the right.
00447     /*
00448     if (Verbose > 0) {
00449         cout << "4/4 ";
00450         cout.flush();
00451     }
00452     */
00453     dist_t *lastDNFLeftP = &dnfLeft[ubbPixels - 1];
00454     m1y = mask1.first + (ubbSize - vigra::Diff2D(1,1));
00455     m2y = mask2.first + (ubbSize - vigra::Diff2D(1,1));
00456     for (uint32 y = 0; y < ubbHeight; y++, --m1y.y, --m2y.y) {
00457         dist_t *dnfColumnP = lastDNFColumnP - (y * ubbWidth);
00458         dist_t *dnfLeftP = lastDNFLeftP - (y * ubbWidth);
00459         //MaskPixel *maskP = lastMaskP - (y * ubbWidth);
00460 
00461         // List of dnfColumnP's on the right that might be the closest features
00462         // to the current dnfColumnP.
00463         std::list<dist_t*> potentialFeatureList;
00464 
00465         Mask1Iter m1x(m1y);
00466         Mask2Iter m2x(m2y);
00467         for (uint32 x = 0; x < ubbWidth; x++, --m1x.x, --m2x.x) {
00468             // First add ourself to the list.
00469             potentialFeatureList.push_back(dnfColumnP);
00470 
00471             // Iterate through list and prune as before.
00472             std::list<dist_t*>::iterator potentialFeature =
00473                     --(potentialFeatureList.end());
00474             dist_t distPotentialFeature = *dnfColumnP;
00475             while (potentialFeature != potentialFeatureList.begin()) {
00476                 // Iterator that points to predecessor.
00477                 std::list<dist_t*>::iterator previousFeature = potentialFeature;
00478                 previousFeature--;
00479 
00480                 // previousFeature is this many columns to the right of (x,y);
00481                 uint32 deltaX = *previousFeature - dnfColumnP;
00482 
00483                 // previousFeature is this far from (x,y);
00484                 dist_t distPreviousFeature = distance(deltaX, *previousFeature);
00485 
00486                 if (distPreviousFeature >= distPotentialFeature) {
00487                     // previousFeature is not a candidate.
00488                     potentialFeatureList.erase(previousFeature);
00489                 } else {
00490                     // previousFeature is a candidate.
00491                     potentialFeature = previousFeature;
00492                     distPotentialFeature = distPreviousFeature;
00493                 }
00494             }
00495 
00496             // The closest feature on the right is potentialFeature.
00497             if (*dnfLeftP > distPotentialFeature) {
00498                 // Recolor maskP.
00499                 //MaskPixel *maskPRight = maskP + (*potentialFeature - dnfColumnP);
00500                 Mask1Iter m1r(m1x);
00501                 int dx = (*potentialFeature - dnfColumnP);
00502                 m1r.x += dx;
00503                 *m1x = *m1r;
00504                 Mask2Iter m2r(m2x);
00505                 m2r.x += dx;
00506                 *m2x = *m2r;
00507             }
00508 
00509             // Move pointers left one column.
00510             dnfLeftP--;
00511             dnfColumnP--;
00512             //maskP--;
00513 
00514             //maxListSize = max(maxListSize, potentialFeatureList.size());
00515         }
00516     }
00517     //if (Verbose > 0) {
00518     //    cout << "max feature list size=" << maxListSize << endl;
00519     //}
00520 
00521     // ensure that only visible pixels are set, apply feat masks to calculated masks
00522     vigra::combineTwoImages(feat1, mask1, mask1,
00523                             vigra::functor::ifThenElse( vigra::functor::Arg1() > vigra::functor::Param(0),
00524                                                         vigra::functor::Arg2(),
00525                                                         vigra::functor::Param(0) ) );
00526     vigra::combineTwoImages(srcIterRange(feat2.first, feat2.first + ubbSize),
00527                             mask2, mask2,
00528                             vigra::functor::ifThenElse( vigra::functor::Arg1() > vigra::functor::Param(0),
00529                                                         vigra::functor::Arg2(),
00530                                                         vigra::functor::Param(0) ) );
00531 
00532 #ifdef DEBUG_NEAREST
00533     vigra::exportImage(vigra::srcIterRange(mask1.first, mask1.first + ubbSize), vigra::ImageExportInfo("nona_4_mask1.tif"));
00534     vigra::exportImage(vigra::srcIterRange(mask2.first, mask2.first + ubbSize), vigra::ImageExportInfo("nona_4_mask2.tif"));
00535     saveDist(dnfColumn, ubbWidth, ubbHeight, "nona_4_dnfColumn.tif");
00536     saveDist(dnfLeft, ubbWidth, ubbHeight, "nona_4_dnfLeft.tif");
00537 #endif
00538     
00539     progress.increase();
00540 
00541     free(dnfColumn);
00542     free(dnfLeft);
00543 
00544     progress.popTask();
00545     return;
00546 }
00547 
00548 } // namespace
00549 
00550 #endif // NEARESTFEATURETRANSFORM_H

Generated on 25 Oct 2014 for Hugintrunk by  doxygen 1.4.7