diff --git a/README.md b/README.md index 1bfc00e..352da44 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ while robust noise and motion estimation maintain broad applicability to many di ### Download -**[Download stable version (0.2.3)](https://github.com/tjof2/pgure-svt/archive/v0.2.3.tar.gz)** +**[Download stable version (0.3.0)](https://github.com/tjof2/pgure-svt/archive/v0.3.0.tar.gz)** PGURE-SVT has been tested on Ubuntu 12.04, 13.04 and 14.04. diff --git a/src/exe-PGURE-SVT.cpp b/src/exe-PGURE-SVT.cpp index 3baf4fa..449c18f 100644 --- a/src/exe-PGURE-SVT.cpp +++ b/src/exe-PGURE-SVT.cpp @@ -1,25 +1,25 @@ /*************************************************************************** - PGURE-SVT Denoising + PGURE-SVT Denoising - Author: Tom Furnival - Email: tjof2@cam.ac.uk + Author: Tom Furnival + Email: tjof2@cam.ac.uk - Copyright (C) 2015-16 Tom Furnival + Copyright (C) 2015-16 Tom Furnival - This program uses Singular Value Thresholding (SVT) [1], combined - with an unbiased risk estimator (PGURE) to denoise a video sequence - of microscopy images [2]. Noise parameters for a mixed Poisson-Gaussian - noise model are automatically estimated during the denoising. + This program uses Singular Value Thresholding (SVT) [1], combined + with an unbiased risk estimator (PGURE) to denoise a video sequence + of microscopy images [2]. Noise parameters for a mixed Poisson-Gaussian + noise model are automatically estimated during the denoising. - References: - [1] "Unbiased Risk Estimates for Singular Value Thresholding and - Spectral Estimators", (2013), Candes, EJ et al. - http://dx.doi.org/10.1109/TSP.2013.2270464 + References: + [1] "Unbiased Risk Estimates for Singular Value Thresholding and + Spectral Estimators", (2013), Candes, EJ et al. + http://dx.doi.org/10.1109/TSP.2013.2270464 - [2] "An Unbiased Risk Estimator for Image Denoising in the Presence - of Mixed Poisson–Gaussian Noise", (2014), Le Montagner, Y et al. - http://dx.doi.org/10.1109/TIP.2014.2300821 + [2] "An Unbiased Risk Estimator for Image Denoising in the Presence + of Mixed Poisson–Gaussian Noise", (2014), Le Montagner, Y et al. + http://dx.doi.org/10.1109/TIP.2014.2300821 This file is part of PGURE-SVT. @@ -66,7 +66,7 @@ namespace libtiff { // Constant-time median filter extern "C" { - #include "medfilter.h" + #include "medfilter.h" } // Own headers @@ -82,95 +82,95 @@ bool strToBool(std::string const& s) {return s != "0";}; // Main program int main(int argc, char** argv) { - // Overall program timer - auto overallstart = std::chrono::steady_clock::now(); + // Overall program timer + auto overallstart = std::chrono::steady_clock::now(); - // Print program header - std::cout< programOptions; - std::ifstream paramFile(argv[1], std::ios::in); - - // Parse the parameter file - ParseParameters(paramFile, programOptions); - - // Check all required parameters are specified - if(programOptions.count("filename") == 0 || programOptions.count("start_image") == 0 || programOptions.count("end_image") == 0) { - std::cout<<"**WARNING** Required parameters not specified"<> tol; - } - - // Block overlap - int Bo = (programOptions.count("patch_overlap") == 1) ? std::stoi(programOptions.at("patch_overlap")) : 1; + std::ifstream paramFile(argv[1], std::ios::in); + + // Parse the parameter file + ParseParameters(paramFile, programOptions); + + // Check all required parameters are specified + if(programOptions.count("filename") == 0 || programOptions.count("start_image") == 0 || programOptions.count("end_image") == 0) { + std::cout<<"**WARNING** Required parameters not specified"<> tol; + } + + // Block overlap + int Bo = (programOptions.count("patch_overlap") == 1) ? std::stoi(programOptions.at("patch_overlap")) : 1; // Noise method // TODO:tjof2 document this option @@ -179,285 +179,285 @@ int main(int argc, char** argv) { // Hot pixel threshold double hotpixelthreshold = (programOptions.count("hot_pixel") == 1) ? std::stoi(programOptions.at("hot_pixel")) : 10; - ///////////////////////////// - // // - // SEQUENCE IMPORT // - // // - ///////////////////////////// - - // Check file exists - std::string infilename = filestem + ".tif"; - if(!std::ifstream(infilename.c_str())) { - std::cout<<"**WARNING** File "<= (startimg-1) && dircount <= (endimg-1)) { - inputsequence.resize(tiffHeight, tiffWidth, imgcount+1); - filteredsequence.resize(tiffHeight, tiffWidth, imgcount+1); - - unsigned short *Buffer = new unsigned short[tiffWidth*tiffHeight]; - unsigned short *FilteredBuffer = new unsigned short[tiffWidth*tiffHeight]; - - for(int tiffRow = 0; tiffRow < tiffHeight; tiffRow++) { - libtiff::TIFFReadScanline(MultiPageTiff, &Buffer[tiffRow*tiffWidth], tiffRow, 0); - } - - arma::Mat TiffSlice( Buffer, tiffHeight, tiffWidth); - inplace_trans(TiffSlice); - inputsequence.slice(imgcount) = arma::conv_to::from(TiffSlice); - - // Apply median filter (constant-time) to the 8-bit image - int memsize = 512 * 1024; // L2 cache size - int filtsize = MedianSize; // Median filter size in pixels - ConstantTimeMedianFilter(Buffer, FilteredBuffer, tiffWidth, tiffHeight, tiffWidth, tiffWidth, filtsize, 1, memsize); - arma::Mat FilteredTiffSlice( FilteredBuffer, tiffHeight, tiffWidth); - inplace_trans(FilteredTiffSlice); - filteredsequence.slice(imgcount) = arma::conv_to::from(FilteredTiffSlice); - imgcount++; - } - dircount++; - } while(libtiff::TIFFReadDirectory(MultiPageTiff)); - libtiff::TIFFClose(MultiPageTiff); - } - // Is number of frames compatible? - if(num_images > (int)inputsequence.n_slices) { - std::cout<<"**WARNING** Sequence only has "<= (startimg-1) && dircount <= (endimg-1)) { + inputsequence.resize(tiffHeight, tiffWidth, imgcount+1); + filteredsequence.resize(tiffHeight, tiffWidth, imgcount+1); + + unsigned short *Buffer = new unsigned short[tiffWidth*tiffHeight]; + unsigned short *FilteredBuffer = new unsigned short[tiffWidth*tiffHeight]; + + for(int tiffRow = 0; tiffRow < tiffHeight; tiffRow++) { + libtiff::TIFFReadScanline(MultiPageTiff, &Buffer[tiffRow*tiffWidth], tiffRow, 0); + } + + arma::Mat TiffSlice( Buffer, tiffHeight, tiffWidth); + inplace_trans(TiffSlice); + inputsequence.slice(imgcount) = arma::conv_to::from(TiffSlice); + + // Apply median filter (constant-time) to the 8-bit image + int memsize = 512 * 1024; // L2 cache size + int filtsize = MedianSize; // Median filter size in pixels + ConstantTimeMedianFilter(Buffer, FilteredBuffer, tiffWidth, tiffHeight, tiffWidth, tiffWidth, filtsize, 1, memsize); + arma::Mat FilteredTiffSlice( FilteredBuffer, tiffHeight, tiffWidth); + inplace_trans(FilteredTiffSlice); + filteredsequence.slice(imgcount) = arma::conv_to::from(FilteredTiffSlice); + imgcount++; + } + dircount++; + } while(libtiff::TIFFReadDirectory(MultiPageTiff)); + libtiff::TIFFClose(MultiPageTiff); + } + // Is number of frames compatible? + if(num_images > (int)inputsequence.n_slices) { + std::cout<<"**WARNING** Sequence only has "<= (num_images - framewindow)) { - u = noisysequence.slices(num_images-2*framewindow-1,num_images-1); - ufilter = filteredsequence.slices(num_images-2*framewindow-1,num_images-1); - } - else { - u = noisysequence.slices(timeiter - framewindow, timeiter + framewindow); - ufilter = filteredsequence.slices(timeiter - framewindow, timeiter + framewindow); - } - - // Basic sequence normalization - double inputmax = u.max(); - u /= inputmax; - ufilter /= ufilter.max(); - - ///////////////////////////// - // // - // NOISE ESTIMATION // - // // - ///////////////////////////// - - // Perform noise estimation - if(pgureOpt) { - NoiseEstimator *noise = new NoiseEstimator; - noise->Estimate(u, - alpha, - mu, - sigma, - 8, - NoiseMethod); - delete noise; - } - - ///////////////////////////// - // // - // MOTION ESTIMATION // - // // - ///////////////////////////// - - MotionEstimator *motion = new MotionEstimator; - motion->Estimate(ufilter, - timeiter, - framewindow, - num_images, - Bs, - MotionP); - arma::icube sequencePatches = motion->GetEstimate(); - delete motion; - - ///////////////////////////// - // // - // PGURE OPTIMIZATION // - // // - ///////////////////////////// - - PGURE *optimizer = new PGURE; - optimizer->Initialize(u, - sequencePatches, - Bs, - Bo, - alpha, - sigma, - mu); - // Determine optimum threshold value (max 1000 evaluations) - if(pgureOpt) { - lambda = (timeiter == 0) ? arma::accu(u)/(Nx*Ny*T) : lambda; - lambda = optimizer->Optimize(tol, lambda, u.max(), 1E3); - v = optimizer->Reconstruct(lambda); - } - else { - v = optimizer->Reconstruct(lambda); - } - delete optimizer; - - ///////////////////////////// - // // - // SEQUENCE RECONSTRUCTION // - // // - ///////////////////////////// - - // Rescale back to original range - v *= inputmax; - - // Place frames back into sequence - if(timeiter < framewindow) { - cleansequence.slice(timeiter) = v.slice(timeiter); - } - else if(timeiter >= (num_images - framewindow)) { - int endseqFrame = timeiter - (num_images - T); - cleansequence.slice(timeiter) = v.slice(endseqFrame); - } - else { - cleansequence.slice(timeiter) = v.slice(framewindow); - } - - // Finish timing - auto endLoopTimer = std::chrono::steady_clock::now(); - auto elapsedLoopTimer = std::chrono::duration_cast(endLoopTimer - startLoopTimer); - - // Output a table row with noise estimates, lambda and timing - std::ostringstream framestring; - framestring << timeiter+1; - - std::cout< outTiff(tiffWidth,tiffHeight,num_images); - outTiff = arma::conv_to>::from(65535*cleansequence); - - // Get the filename - std::string outfilename = filestem + "-CLEANED.tif"; - - // Set the output file headers - libtiff::TIFF *MultiPageTiffOut = libtiff::TIFFOpen(outfilename.c_str(), "w"); - libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_IMAGEWIDTH, tiffWidth); - libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_IMAGELENGTH, tiffHeight); - libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_BITSPERSAMPLE, 16); - - // Write the file - if(!MultiPageTiffOut) { - std::cout<<"**WARNING** File "<= (num_images - framewindow)) { + u = noisysequence.slices(num_images-2*framewindow-1,num_images-1); + ufilter = filteredsequence.slices(num_images-2*framewindow-1,num_images-1); + } + else { + u = noisysequence.slices(timeiter - framewindow, timeiter + framewindow); + ufilter = filteredsequence.slices(timeiter - framewindow, timeiter + framewindow); + } + + // Basic sequence normalization + double inputmax = u.max(); + u /= inputmax; + ufilter /= ufilter.max(); + + ///////////////////////////// + // // + // NOISE ESTIMATION // + // // + ///////////////////////////// + + // Perform noise estimation + if(pgureOpt) { + NoiseEstimator *noise = new NoiseEstimator; + noise->Estimate(u, + alpha, + mu, + sigma, + 8, + NoiseMethod); + delete noise; + } + + ///////////////////////////// + // // + // MOTION ESTIMATION // + // // + ///////////////////////////// + + MotionEstimator *motion = new MotionEstimator; + motion->Estimate(ufilter, + timeiter, + framewindow, + num_images, + Bs, + MotionP); + arma::icube sequencePatches = motion->GetEstimate(); + delete motion; + + ///////////////////////////// + // // + // PGURE OPTIMIZATION // + // // + ///////////////////////////// + + PGURE *optimizer = new PGURE; + optimizer->Initialize(u, + sequencePatches, + Bs, + Bo, + alpha, + sigma, + mu); + // Determine optimum threshold value (max 1000 evaluations) + if(pgureOpt) { + lambda = (timeiter == 0) ? arma::accu(u)/(Nx*Ny*T) : lambda; + lambda = optimizer->Optimize(tol, lambda, u.max(), 1E3); + v = optimizer->Reconstruct(lambda); + } + else { + v = optimizer->Reconstruct(lambda); + } + delete optimizer; + + ///////////////////////////// + // // + // SEQUENCE RECONSTRUCTION // + // // + ///////////////////////////// + + // Rescale back to original range + v *= inputmax; + + // Place frames back into sequence + if(timeiter < framewindow) { + cleansequence.slice(timeiter) = v.slice(timeiter); + } + else if(timeiter >= (num_images - framewindow)) { + int endseqFrame = timeiter - (num_images - T); + cleansequence.slice(timeiter) = v.slice(endseqFrame); + } + else { + cleansequence.slice(timeiter) = v.slice(framewindow); + } + + // Finish timing + auto endLoopTimer = std::chrono::steady_clock::now(); + auto elapsedLoopTimer = std::chrono::duration_cast(endLoopTimer - startLoopTimer); + + // Output a table row with noise estimates, lambda and timing + std::ostringstream framestring; + framestring << timeiter+1; + + std::cout< outTiff(tiffWidth,tiffHeight,num_images); + outTiff = arma::conv_to>::from(65535*cleansequence); + + // Get the filename + std::string outfilename = filestem + "-CLEANED.tif"; + + // Set the output file headers + libtiff::TIFF *MultiPageTiffOut = libtiff::TIFFOpen(outfilename.c_str(), "w"); + libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_IMAGEWIDTH, tiffWidth); + libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_IMAGELENGTH, tiffHeight); + libtiff::TIFFSetField(MultiPageTiffOut, TIFFTAG_BITSPERSAMPLE, 16); + + // Write the file + if(!MultiPageTiffOut) { + std::cout<<"**WARNING** File "< outSlice = outTiff.slice(tOut); - inplace_trans(outSlice); - unsigned short *OutBuffer = outSlice.memptr(); - libtiff::TIFFWriteScanline(MultiPageTiffOut, &OutBuffer[tiffRow*tiffWidth], tiffRow, 0); - } - libtiff::TIFFWriteDirectory(MultiPageTiffOut); - } - libtiff::TIFFClose(MultiPageTiffOut); - - ///////////////////////////// - // // - // REPORT RESULT // - // // - ///////////////////////////// - - // Overall program timer - auto overallend = std::chrono::steady_clock::now(); - auto elapsed = std::chrono::duration_cast(overallend - overallstart); - std::cout<<"Total time: "< outSlice = outTiff.slice(tOut); + inplace_trans(outSlice); + unsigned short *OutBuffer = outSlice.memptr(); + libtiff::TIFFWriteScanline(MultiPageTiffOut, &OutBuffer[tiffRow*tiffWidth], tiffRow, 0); + } + libtiff::TIFFWriteDirectory(MultiPageTiffOut); + } + libtiff::TIFFClose(MultiPageTiffOut); + + ///////////////////////////// + // // + // REPORT RESULT // + // // + ///////////////////////////// + + // Overall program timer + auto overallend = std::chrono::steady_clock::now(); + auto elapsed = std::chrono::duration_cast(overallend - overallstart); + std::cout<<"Total time: "<= 0.) ? userLambda : 0.;