[Insight-users] Color Deconvolution in ITK for H&E stained images
Dan Mueller
dan.muel at gmail.com
Tue Dec 11 18:18:18 EST 2012
Hi Humayun,
During a previous project I needed to segment some HE histology
images. I ended up being able to simply use the greyscale image, but
as part of the process I briefly investigated so-called "color
deconvolution" (if I recall correctly, I only spent about a day or two
investigating).
I managed to dig up my old code (see below). The difficult part of
automatically computing the OpticalDensityMatrix is not implemented,
but the code as-is might provide some inspiration. WARNING: I do not
guarantee this code in way whatsoever.
The source comments have the following URLs which it seems I found
helpful at the time:
http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
http://svn.openmicroscopy.org.uk/svn/ome/trunk/src/matlab/OME/Filters/colour_deconvolution/colour_deconvolution.c
https://bitbucket.org/David_Clark/color-deconvolution/src/1b8b24e6f59ec1c54c40d90a88a72b02bb777e7b/Reference/Quantification%20of%20Histochemical%20Staining%20by%20Color%20Deconvolution%20-%20%20Ruifrok%2C%20A.pdf
HTH
Cheers, Dan
== CMakeLists.txt ==
# This project is designed to be built outside the Insight source tree.
PROJECT(ColorDeconvolution)
CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
# Find ITK
FIND_PACKAGE(ITK REQUIRED)
INCLUDE(${ITK_USE_FILE})
# Find Boost
SET(Boost_USE_STATIC_LIBS ON)
FIND_PACKAGE(Boost 1.51 REQUIRED)
INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIR})
LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
# The header files
SET(HEADERS
itkColorDeconvolutionImageFilter.h
)
# The implementation files
SET(SOURCES
ColorDeconvolution.cxx
)
# Set include directories
INCLUDE_DIRECTORIES(
${INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}
${SOURCE_PATH}
${BOOST_INCLUDE_DIR}
)
# Set link directories
LINK_DIRECTORIES(
${LINK_DIRECTORIES}
${BOOST_LINK_DIR}
)
# Main library
ADD_EXECUTABLE(
ColorDeconvolution
${HEADERS}
${SOURCES}
)
TARGET_LINK_LIBRARIES(
ColorDeconvolution
${ITK_LIBRARIES}
${Boost_LIBRARIES}
)
== ColorDeconvolution.cxx ==
//
// Filename: ColorDeconvolution.cxx
//
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#define _SCL_SECURE_NO_WARNINGS
#define NO_TESTING
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdio>
#include "boost/program_options.hpp"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkColorDeconvolutionImageFilter.h"
#include "itkUnaryFunctorImageFilter.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
// Helpful typedefs
const unsigned int Dimensions = 2;
const unsigned int Channels = 3;
typedef unsigned char ScalarPixelType;
typedef itk::Vector<ScalarPixelType, Channels> VectorPixelType;
typedef itk::Image<ScalarPixelType, Dimensions> ScalarImageType;
typedef itk::Image<VectorPixelType, Dimensions> VectorImageType;
typedef itk::ImageFileReader<ScalarImageType> ScalarReaderType;
typedef itk::ImageFileWriter<ScalarImageType> ScalarWriterType;
typedef itk::ImageFileReader<VectorImageType> VectorReaderType;
typedef itk::ImageFileWriter<VectorImageType> VectorWriterType;
typedef itk::ColorDeconvolutionImageFilter<VectorImageType>
ColorDeconvolutionType;
typedef ColorDeconvolutionType::OpticalDensityMatrixType
OpticalDensityMatrixType;
typedef itk::Functor::ColorDeconvolutionLookup<ScalarPixelType,
VectorPixelType> LookupFunctorType;
typedef itk::UnaryFunctorImageFilter<ScalarImageType, VectorImageType,
LookupFunctorType> LookupImageFilterType;
typedef itk::VectorIndexSelectionCastImageFilter<VectorImageType,
ScalarImageType> VectorIndexSelectType;
// Function prototypes
std::string ConstructPostfixedFilename(std::string name, std::string
postfix, std::string ext);
int main(int argc, char * argv []) {
try {
// Parse command line arguments
std::string inputFilename;
std::string outputFilename;
boost::program_options::options_description options(
"Usage: ColorDeconvolution [options] input-file
output-file:\n\noptions");
boost::program_options::variables_map inputMap;
{
using namespace boost::program_options;
options.add_options()
("help", "print this help message")
("input-file,i" ,
value<std::string>(&inputFilename)->required(),
"input image filename")
("output-file,o",
value<std::string>(&outputFilename)->required(),
"output csv filename")
;
positional_options_description positionalOptions;
positionalOptions.add("input-file", 1);
positionalOptions.add("output-file", 1);
try {
store(command_line_parser(argc, argv).options(options).
positional(positionalOptions).run(),inputMap);
notify(inputMap);
} catch (...) {
std::cout << "Error: invalid options specified.\n\n"
<< options << "\n";
return EXIT_FAILURE;
}
}
// Read input image
std::cout << "Reading input image: " << inputFilename << std::endl;
VectorReaderType::Pointer reader = VectorReaderType::New();
reader->SetFileName(inputFilename);
reader->Update();
// Apply color deconvolution
ColorDeconvolutionType::Pointer deconvolution =
ColorDeconvolutionType::New();
deconvolution->SetInput(reader->GetOutput());
OpticalDensityMatrixType odm;
// Haematoxylin and Eosin separation
odm[0][0] = 0.65;
odm[0][1] = 0.70;
odm[0][2] = 0.29;
odm[1][0] = 0.07;
odm[1][1] = 0.99;
odm[1][2] = 0.11;
odm[2][0] = 0.27;
odm[2][1] = 0.57;
odm[2][2] = 0.78;
odm[0][0] = 0.18;
odm[0][1] = 0.20;
odm[0][2] = 0.08;
odm[1][0] = 0.01;
odm[1][1] = 0.13;
odm[1][2] = 0.01;
odm[2][0] = 0.10;
odm[2][1] = 0.21;
odm[2][2] = 0.29;
deconvolution->SetOpticalDensityMatrix(odm);
deconvolution->Update();
VectorImageType::Pointer deconvolved = deconvolution->GetOutput();
deconvolved->DisconnectPipeline();
// Create color images from deconvolved channel
for (unsigned int i=0; i<Channels; i++) {
// Select channel
VectorIndexSelectType::Pointer select =
VectorIndexSelectType::New();
select->SetInput(deconvolved);
select->SetIndex(i);
select->Update();
// Apply color LUT
LookupImageFilterType::Pointer lut = LookupImageFilterType::New();
lut->SetInput(select->GetOutput());
LookupFunctorType::OpticalDensityVectorType odv;
for (unsigned int c=0; c<Channels; c++) {
odv[c] = deconvolution->GetOpticalDensityMatrix()[i][c];
}
lut->GetFunctor().SetOpticalDensityVector(odv);
lut->Update();
// Write
char channelName[128];
std::sprintf(channelName, "-channel%d", i);
std::string channelFilename = ConstructPostfixedFilename(
outputFilename, channelName, ".png"
);
std::cout << "Writing channel image: " << channelFilename
<< std::endl;
VectorWriterType::Pointer writer = VectorWriterType::New();
writer->SetFileName(channelFilename);
writer->SetInput(lut->GetOutput());
writer->Update();
}
// Write output
std::cout << "Writing output image: " << outputFilename << std::endl;
VectorWriterType::Pointer writer = VectorWriterType::New();
writer->SetFileName(outputFilename);
writer->SetInput(deconvolved);
writer->Update();
return EXIT_SUCCESS;
}
catch (itk::ExceptionObject &e) {
std::cerr << e << std::endl;
return EXIT_FAILURE;
}
}
// Construct filename with the given postfix inserted between the
filename and given extension.
// "file" could be a fully qualified or relative path.
// "ext" must include the dot. If not specified, the given filename
extension is used.
// Example1: file="C:/A/B/C.ext", postfix="-D", ext=".png",
result="C:/A/B/C-D.png"
// Example2: file="C:/A/B/C.ext", postfix="-D", ext="",
result="C:/A/B/C-D.ext"
std::string ConstructPostfixedFilename(std::string file, std::string
postfix, std::string ext) {
std::string path = itksys::SystemTools::GetFilenamePath(file);
std::string fileWithoutExtension =
itksys::SystemTools::GetFilenameWithoutLastExtension(file);
std::string fileExtension =
itksys::SystemTools::GetFilenameLastExtension(file);
std::string result = path + "/" + fileWithoutExtension + postfix;
result += (ext.length() > 0) ? ext : fileExtension;
return result;
}
== itkColorDeconvolutionImageFilter.h ==
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef __itkColorDeconvolutionImageFilter_h
#define __itkColorDeconvolutionImageFilter_h
#include "vnl/vnl_math.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_inverse.h"
#include "itkVector.h"
#include "itkUnaryFunctorImageFilter.h"
namespace itk
{
/** Functor for color deconvolution. */
namespace Functor {
template<class TPixel>
class ColorDeconvolution
{
public:
typedef TPixel PixelType;
typedef typename TPixel::ValueType ValueType;
typedef typename NumericTraits<ValueType>::RealType RealType;
typedef vnl_matrix_fixed<RealType, 3, 3> OpticalDensityMatrixType;
typedef vnl_matrix_fixed<RealType, 3, 3> ColorDeconvolutionMatrixType;
ColorDeconvolution() {}
~ColorDeconvolution() {}
virtual const char *GetNameOfClass() const
{
return "Functor::ColorDeconvolution";
}
bool operator!=( const ColorDeconvolution & ) const
{
return false;
}
bool operator==( const ColorDeconvolution & other ) const
{
return !(*this != other);
}
typename OpticalDensityMatrixType GetOpticalDensityMatrix() const
{
return m_OpticalDensityMatrix;
}
void SetOpticalDensityMatrix(typename OpticalDensityMatrixType odm)
{
odm.normalize_rows();
// If a third color is not specified, compute the compliment
if ((odm[2][0] == 0.0) && (odm[2][1] == 0.0) && (odm[2][2] == 0.0))
{
odm[2][0] = ((odm[0][0]*odm[0][0] + odm[1][0]*odm[1][0]) > 1.0) ?
0.0 : sqrt(1.0-(odm[0][0]*odm[0][0])-(odm[1][0]*odm[1][0]));
odm[2][1] = ((odm[0][1]*odm[0][1] + odm[1][1]*odm[1][1]) > 1.0) ?
0.0 : sqrt(1.0-(odm[0][1]*odm[0][1])-(odm[1][1]*odm[1][1]));
odm[2][2] = ((odm[0][2]*odm[0][2] + odm[1][2]*odm[1][2]) > 1.0) ?
0.0 : sqrt(1.0-(odm[0][2]*odm[0][2])-(odm[1][2]*odm[1][2]));
}
// Renormalize, for any changes in last row
odm.normalize_rows();
// Check for zero values
for (unsigned int i=0; i<3; i++)
{
for (unsigned int j=0; j<3; j++)
{
odm[i][j] = (odm[i][j] == 0.0) ? 0.0001 : odm[i][j];
}
}
m_OpticalDensityMatrix = odm;
// Invert to create color deconvolution matrix
m_ColorDeconvolutionMatrix = vnl_inverse(odm).transpose();
}
inline TPixel operator()( const TPixel & input ) const
{
TPixel output;
Vector<RealType, PixelType::Dimension> logValues;
RealType log255 = log(255.0);
RealType scaled = 0.0;
for (unsigned int i=0; i<3; i++)
{
logValues[i] =
-((255.0*log((static_cast<RealType>(input[i]))/255.0))/log255);
}
for (unsigned int i=0; i<3; i++)
{
scaled = 0.0;
for (unsigned int j=0; j<3; j++)
{
scaled += logValues[j] * m_ColorDeconvolutionMatrix[i][j];
}
double outputReal = exp(-(scaled-255.0) * log255 / 255.0);
output[i] = (outputReal > NumericTraits<ValueType>::max()) ?
NumericTraits<ValueType>::max() : static_cast<ValueType>(outputReal);
}
return output;
}
private:
OpticalDensityMatrixType m_OpticalDensityMatrix;
ColorDeconvolutionMatrixType m_ColorDeconvolutionMatrix;
};
/** Functor for applying RGB LUT to single channel image. */
template<class TInput, class TOutput>
class ColorDeconvolutionLookup
{
public:
typedef TInput InputPixelType;
typedef TOutput OutputPixelType;
typedef typename OutputPixelType::ValueType OutputValueType;
typedef typename NumericTraits<InputPixelType>::RealType RealType;
typedef FixedArray<RealType, 3> OpticalDensityVectorType;
ColorDeconvolutionLookup() {}
~ColorDeconvolutionLookup() {}
virtual const char *GetNameOfClass() const
{
return "Functor::ColorDeconvolutionLookup";
}
bool operator!=( const ColorDeconvolutionLookup & ) const
{
return false;
}
bool operator==( const ColorDeconvolutionLookup & other ) const
{
return !(*this != other);
}
typename OpticalDensityVectorType GetOpticalDensityVector() const
{
return m_OpticalDensityVector;
}
void SetOpticalDensityVector(const typename OpticalDensityVectorType odv)
{
m_OpticalDensityVector = odv;
}
inline TOutput operator()( const TInput & input ) const
{
TOutput output;
OutputValueType temp;
for (unsigned int i=0; i<OutputPixelType::Dimension; i++)
{
temp = 255.0 - ((255.0-static_cast<RealType>(input)) *
m_OpticalDensityVector[i]);
output[i] = static_cast<OutputValueType>(temp);
}
return output;
}
private:
OpticalDensityVectorType m_OpticalDensityVector;
};
} // end namespace Functor
/** \class ColorDeconvolutionImageFilter
* \brief "Deconvolves" color spaces.
*
* \see http://www.dentistry.bham.ac.uk/landinig/software/cdeconv/cdeconv.html
* \see http://svn.openmicroscopy.org.uk/svn/ome/trunk/src/matlab/OME/Filters/colour_deconvolution/colour_deconvolution.c
* \see https://lcib.rutgers.edu/~james/quantcolor.pdf
* \see https://bitbucket.org/David_Clark/color-deconvolution/src/1b8b24e6f59ec1c54c40d90a88a72b02bb777e7b/Reference/Quantification%20of%20Histochemical%20Staining%20by%20Color%20Deconvolution%20-%20%20Ruifrok%2C%20A.pdf
*
* \ingroup IntensityImageFilters Multithreaded
*/
template <class TImage>
class ITK_EXPORT ColorDeconvolutionImageFilter :
public UnaryFunctorImageFilter<TImage,TImage,
Functor::ColorDeconvolution<typename TImage::PixelType> >
{
public:
/** Standard class typedefs. */
typedef ColorDeconvolutionImageFilter Self;
typedef UnaryFunctorImageFilter<TImage, TImage,
Functor::ColorDeconvolution<typename TImage::PixelType> > Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Runtime information support. */
itkTypeMacro(ColorDeconvolutionImageFilter,
UnaryFunctorImageFilter);
typedef typename TImage::PixelType PixelType;
typedef typename NumericTraits<PixelType>::ValueType PixelValueType;
typedef typename Functor::ColorDeconvolution<PixelType> FunctorType;
typedef typename FunctorType::OpticalDensityMatrixType
OpticalDensityMatrixType;
/** Get/set OpticalDensityMatrix. */
virtual OpticalDensityMatrixType GetOpticalDensityMatrix() const
{
return this->GetFunctor().GetOpticalDensityMatrix();
}
virtual void SetOpticalDensityMatrix(const OpticalDensityMatrixType odm)
{
if( odm != this->GetFunctor().GetOpticalDensityMatrix() )
{
this->GetFunctor().SetOpticalDensityMatrix( odm );
this->Modified();
}
}
#ifdef ITK_USE_CONCEPT_CHECKING
/** Begin concept checking */
// TODO:
/** End concept checking */
#endif
protected:
ColorDeconvolutionImageFilter() {}
virtual ~ColorDeconvolutionImageFilter() {}
private:
ColorDeconvolutionImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
} // end namespace itk
#endif
On 11 December 2012 12:21, Humayun Irshad <humayun.irshad at gmail.com> wrote:
>
> Dear Insight Users,
>
> I am searching a filer in ITK for color deconvolution. I found following two
> article for deconvolution:
>
> Noise Simulation (http://www.insight-journal.org/browse/publication/721)
> FFT based Convolution
> (http://www.insight-journal.org/download/viewpdf/717/4)
>
> Both of these article are related with image debluring. I am working on
> stain separation for histological images, and I was wondering if there are
> by any chance already implemented routines for color deconvolution? Can you
> help me regarding it.
>
> --
> --
> Best Regards,
>
> HUMAYUN IRSHAD
> PhD student, University of Joseph Fourier, Grenoble, France
> Research Engineer,
> IPAL – Image & Pervasive Access Lab, Singapore
> UMI CNRS (I2R/A*STAR, NUS, UJF, UPMC, IT)
>
> Institute for Infocomm Research (I2R)
> 1 Fusionopolis Way
> #10-19 Connexis South Tower
> Singapore 138632
> Lab Tel:+65-64082601
> Cell: +65-86573510
> +60-1128642416
>
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list