PreprocessingSentinel.cxx

/*
 * Copyright (C) 2018-2019, Centre National d'Etudes Spatiales (CNES)
 * All rights reserved
 *
 * This file is part of Weighted Average Synthesis Processor (WASP)
 *
 * Authors:
 * - Peter KETTIG <peter.kettig@cnes.fr>
 *
 * 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 3 of the License, or (at
 * your option) any later version.
 *
 * See the LICENSE.md file for more details.
 *
 * 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, see <https://www.gnu.org/licenses/>.
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */

#include "PreprocessingSentinel.h"
#include "PreprocessingAdapter.h"
#include "DirectionalCorrection.h"
#include "MetadataHelperFactory.h"

#include "otbWrapperApplication.h"

using namespace ts::preprocessing;

const char * PreprocessingSentinel::GetNameOfClass() const
{
  return "PreprocessingSentinel";
}

PreprocessingAdapter::FloatImageType::Pointer PreprocessingSentinel::getCloudMask(const std::string &filename, const unsigned char &bit){
	return getFloatMask(filename, bit);
}

PreprocessingAdapter::FloatImageType::Pointer PreprocessingSentinel::getAotMask(const std::string &filename, const unsigned char &band){

	m_aotReader = PreprocessingAdapter::FloatVectorImageReaderType::New();
	m_aotReader->SetFileName(filename);
	return m_aotExtractor.ExtractImgResampledBand(m_aotReader->GetOutput(), band, Interpolator_Linear);
}

PreprocessingAdapter::FloatImageType::Pointer PreprocessingSentinel::getWaterMask(const std::string &filename, const unsigned char &bit ){
	return getFloatMask(filename, bit);
}

PreprocessingAdapter::FloatImageType::Pointer PreprocessingSentinel::getSnowMask(const std::string &filename, const unsigned char &bit){
	return getFloatMask(filename, bit);
}

std::vector<PreprocessingAdapter::ShortVectorImageType::Pointer> PreprocessingSentinel::getCorrectedRasters(
		const std::string &filename,
		FloatImageType::Pointer cloudImage, FloatImageType::Pointer watImage,
		FloatImageType::Pointer snowImage){

	std::vector<PreprocessingAdapter::ShortVectorImageType::Pointer> outputRasters;

	auto factory = ts::MetadataHelperFactory::New();
	auto pHelper = factory->GetMetadataHelper(filename);

	if(m_DirCorMode == ROY && m_scatteringCoeffs.empty()){
		itkExceptionMacro("Need to set scattering coefficients before running correction for S2.");
	}
	if(m_DirCorMode == LUT && m_DirCorLutXmlFileName.empty()){
			itkExceptionMacro("Need to set lut xml file before running LUT directional correction correction for S2.");
		}

	// Compute NDVI for Primary resolution
	PreprocessingAdapter::FloatImageType::Pointer ndviImg;
	if (m_DirCorMode != LUT){
		m_computeNdvi.DoInit(filename);
		ndviImg = m_computeNdvi.DoExecute();
	}

	size_t totalNRes = pHelper->getResolutions().getNumberOfResolutions();

	//For all possible resolutions, do the following
	for(size_t resolution = 0; resolution < totalNRes; resolution++){
		CreateS2AnglesRaster createAngles;
		createAngles.DoInit(resolution, filename);
		PreprocessingAdapter::FloatVectorImageType::Pointer anglesImg = createAngles.DoExecute();
		anglesImg->Update();
		m_createAngles.push_back(createAngles);

		//If Resolution is not principal Resolution, then resize the additional images
		std::cout << "Current resolution: " <<
				pHelper->getResolutions().getResolutionVector()[resolution].getBands()[0].getResolution() << std::endl;
		PreprocessingAdapter::FloatImageType::Pointer ndviImgToUse;
		PreprocessingAdapter::FloatImageType::Pointer cldImgToUse;
		PreprocessingAdapter::FloatImageType::Pointer watImgToUse;
		PreprocessingAdapter::FloatImageType::Pointer snowImgToUse;
		if(cloudImage.GetPointer()->GetSignedSpacing()[0] !=
				pHelper->getResolutions().getResolutionVector()[resolution].getBands()[0].getResolution()){
			if (m_DirCorMode != LUT){
				ndviImgToUse =	m_Resampler.getResampler(ndviImg.GetPointer(), 0.5f)->GetOutput();
			}
			cldImgToUse = m_Resampler.getResampler(cloudImage.GetPointer(), 0.5f)->GetOutput();
			watImgToUse = m_Resampler.getResampler(watImage.GetPointer(), 0.5f)->GetOutput();
			snowImgToUse = m_Resampler.getResampler(snowImage.GetPointer(), 0.5f)->GetOutput();
		} else {
			if (m_DirCorMode != LUT){
				ndviImgToUse =	ndviImg;
			}
			cldImgToUse = cloudImage;
			watImgToUse = watImage;
			snowImgToUse = snowImage;
		}
		ts::DirectionalCorrection dirCorr;
		ts::DirectionalCorrectionLUT dirCorrLut;
		switch (m_DirCorMode) {
			case ROY:
				dirCorr.Init(resolution, filename, m_scatteringCoeffs[resolution],
						cldImgToUse, watImgToUse, snowImgToUse, anglesImg, ndviImgToUse);
				dirCorr.DoExecute();
				m_dirCorr.push_back(dirCorr);
				outputRasters.push_back(dirCorr.GetCorrectedImg().GetPointer());
				break;
			case LUT:
				dirCorrLut.Init(resolution, filename, m_DirCorLutXmlFileName,
						cldImgToUse, watImgToUse, snowImgToUse, anglesImg);
				dirCorrLut.DoExecute();
				m_dirCorrLut.push_back(dirCorrLut);
				outputRasters.push_back(dirCorrLut.GetCorrectedImg().GetPointer());
				break;
			default:
				break;
		}

	}
	return outputRasters;
}