qualicheck.cpp

#include "qualicheck.h"
#include "ui_qualicheck.h"
#include <algorithm>
#include <isis/data/io_application.hpp>
#include <isis/adapter/qt5/simpleimageview.hpp>
#include <isis/data/io_factory.hpp>
#include <isis/math/histogram.hpp>
#include <gsl/gsl_linalg.h>
#include <numeric>

using namespace isis;

// std::array<size_t,256> mk_histogram(data::TypedChunk<uint8_t> image){
// 	std::array<size_t,256> ret;
// 	ret.fill(0);
// 	for(uint8_t v:image)
// 		ret[v]++;
// 	return ret;
// }
// 
// std::array<double,256> mk_normalized_histogram(data::TypedChunk<uint8_t> image){
// 	std::array<double,256> ret;
// 	const std::array<size_t,256> hist=mk_histogram(image);
// 	const uint64_t sum=std::accumulate(hist.begin(),hist.end(),0);
// 
// 	
// 	std::transform(hist.begin(),hist.end(),ret.begin(),[sum](size_t v){return double(v)/sum;});
// 
// 	return ret;
// }

uint8_t otsu_thres(data::TypedChunk<uint8_t> image){
	std::array<double,256> p=math::normalized_histogram(image);

	//M: omega = cumsum(p);
	std::array<double,256> omega;
	std::partial_sum(p.begin(),p.end(),omega.begin()); 
	
	//M: mu = cumsum(p .* (1:num_bins)');
	std::array<double,256> range,prod,mu;
	std::iota(range.begin(),range.end(),1); 
	prod = p * range;
	std::partial_sum(prod.begin(),prod.end(),mu.begin()); 
	double mu_t = mu.back();
	
	//M:sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));
	std::array<double,256> one;one.fill(1);
	const auto v1 = (omega*mu_t - mu) * (omega*mu_t - mu);
	const auto v2 = omega * (one - omega);
	std::array<double,256> sigma_b_squared = v1 / v2;

	double max = *std::max_element(sigma_b_squared.begin(),sigma_b_squared.end());
	
	if(std::isfinite(max)){
		//M: idx = mean(find(sigma_b_squared == maxval));
		const size_t idx1 = std::distance(sigma_b_squared.begin(), std::find(sigma_b_squared.begin(),sigma_b_squared.end(),max));
		const size_t idx2 = std::distance(std::find(sigma_b_squared.rbegin(),sigma_b_squared.rend(),max), sigma_b_squared.rend());
		const size_t idx = idx1+ (idx2-idx1)/2;
		return idx;
	} else
		return 0;
}

data::TypedChunk<bool> binarise(data::TypedChunk<uint8_t> ch){
	uint8_t thres=otsu_thres(ch);
	data::scaling_pair scale(util::Value<int>(1),util::Value<int>(-thres));
	return ch.copyByID(data::ValueArray<bool>::staticID(),scale);
}

std::list<std::pair<size_t,size_t>> traceboundary(data::TypedChunk<bool> ch, std::pair<size_t,size_t> start){ 
	//http://www.imageprocessingplace.com/downloads_V3/root_downloads/tutorials/contour_tracing_Abeer_George_Ghuneim/square.html
	std::list<std::pair<size_t,size_t>> ret;
	struct tracer:std::pair<size_t,size_t>{
		tracer(const std::pair<size_t,size_t> &p):std::pair<size_t,size_t>(p),direction(north){}
		enum D{north=0,east,south,west}direction;
		void move(){
			switch(direction){
				case north:second++;break;
				case east:first++;break;
				case south:second--;break;
				case west:first--;break;
			}
		}
		void go_left(){
			direction = D((uint(direction)+1)%4);
			move();
		}
		void go_right(){
			direction = D((uint(direction)-1)%4);
			move();
		}
	}p(start);
	
	do{
		if(ch.voxel<bool>(p.first,p.second)){//turn left if black
			ret.push_back(p);
			p.go_left();
		} else { // right otherwise
			p.go_right();
		}
	}while(p!=start);
	return ret;
}

std::tuple<double,double,double> solve_gsl(std::list<std::pair<size_t,size_t>> boundary){

	gsl_matrix *m = gsl_matrix_alloc(boundary.size(), 3);

	gsl_vector *b = gsl_vector_alloc(boundary.size());

	gsl_vector *x = gsl_vector_alloc (3);
	
	size_t i=0;
	for(auto p:boundary){ 
		p.first++;
		p.second++;
		//M:[x y ones(length(x),1)]
		gsl_matrix_set(m,i,0,p.first);
		gsl_matrix_set(m,i,1,p.second);
		gsl_matrix_set(m,i,2,1);
		
		//M:[-(x.^2+y.^2)]
		gsl_vector_set(b,i,-double(p.first*p.first+p.second*p.second));
		i++;
	}

	gsl_vector *residual = gsl_vector_alloc (boundary.size());
	gsl_vector *tau = gsl_vector_alloc (3);
	gsl_linalg_QR_decomp (m, tau); // 
	gsl_linalg_QR_lssolve(m, tau, b, x, residual);

	printf ("x = \n");
	gsl_vector_fprintf (stdout, x, "%g");

	return std::make_tuple(gsl_vector_get(x,0),gsl_vector_get(x,1),gsl_vector_get(x,2));
}

template<typename T1, typename T2> data::TypedChunk<T1>& operator +=(data::TypedChunk<T1> &ch1, data::TypedChunk<T2> ch2){
	std::transform(
		ch1.begin(),ch1.end(),
		ch2.begin(),
		ch1.begin(),
		[](const T1 &v,const T2 &vodd)->T1{return v+vodd;}
	);
	return ch1;
}
template<typename T1, typename T2> data::TypedChunk<T1>& operator +=(data::TypedChunk<T1> &ch1, T2 scalar){
	for(T1 &v:ch1)
		v+=scalar;
	return ch1;
}

template<typename T1, typename T2> data::TypedChunk<decltype(T1()+T2())> operator +(data::TypedChunk<T1> ch1, data::TypedChunk<T2> ch2){
	data::MemChunk<decltype(T1()+T2())> ret(ch1);
	return ret+=ch2;
}

template<typename T1, typename T2> data::TypedChunk<decltype(T1()*T2())> operator *(data::TypedChunk<T1> ch1, data::TypedChunk<T2> ch2){
	const auto size = ch1.getSizeAsVector();
	data::MemChunk<decltype(T1()*T2())> ret(size[0],size[1],size[2],size[3]);
	static_cast<util::PropertyMap&>(ret)=ch1;
	std::transform(
		ch1.begin(),ch1.end(),
		ch2.begin(),
		ret.begin(),
		[](const T1 &v,const T2 &vodd)->decltype(T1()*T2()){return v*vodd;}
	);
	return ret;
}
template<typename T1, typename T2> data::TypedChunk<decltype(T1()-T2())> operator -(data::TypedChunk<T1> ch1, data::TypedChunk<T2> ch2){
	const auto size = ch1.getSizeAsVector();
	data::MemChunk<decltype(T1()-T2())> ret(size[0],size[1],size[2],size[3]);
	static_cast<util::PropertyMap&>(ret)=ch1;
	std::transform(
		ch1.begin(),ch1.end(),
		ch2.begin(),
		ret.begin(),
		[](const T1 &v1,const T2 &v2)->decltype(T1()-T2()){return v1-v2;}
	);
	return ret;
}

template<typename T1, typename T2> data::TypedChunk<decltype(T1()*T2())> operator *(data::TypedChunk<T1> ch1, T2 scalar){
	data::MemChunk<decltype(T1()*T2())> ret(ch1);
	for(auto &v:ret)
		v*=scalar;
	return ret;
}

qualicheck::qualicheck(std::list<data::Image> images, QWidget *parent) :  QDialog(parent),  ui(new Ui::Dialog)
{
    ui->setupUi(this);

	data::Image src=images.front();
	
	ui->image_frame->layout()->addWidget(new isis::qt5::SimpleImageView(src,"source"));
	
	src.spliceDownTo(data::sliceDim);
	size_t slices=src.getNrOfSlices();
	
	const double R = src.getNrOfColumns()==128?30:15;
	const std::array<size_t, 4> slize_size=src.getChunkAt(0).getSizeAsVector();
	
	for(size_t s=0;s<slices;s++){
		data::Chunk ch=binarise(src.getChunk(0,0,s,0));
		size_t col = 79, row = 0;
		for(;ch.voxel<bool>(col,row,0)==false && row<slize_size[0];row++){}
		int connectivity = 8;
		size_t num_points   = 180;
		auto contour=traceboundary(ch,{col,row});
		std::cout << contour << std::endl;
		
		double a,b,c;
		std::tie(a, b, c)=solve_gsl(contour);
		
		//M: calculate the location of the center and the radius
		size_t xc = std::round(-a/2),yc = std::round(-b/2);
		double radius  =  sqrt((xc*xc+yc*yc)-c);

		//M: ROI Koordinaten
		size_t X1 = xc - std::floor(R/2);
		size_t Y1 = yc - std::floor(R/2); 
		
		size_t S0=0,St=0,Stt=0;
		std::vector<double> roi;

		data::MemChunk<int> Iodd(slize_size[0],slize_size[1],1,1,true),Ieven(slize_size[0],slize_size[1],1,1,true);
		data::MemChunk<double> Syt(slize_size[0],slize_size[1]),Syy(slize_size[0],slize_size[1]);
		std::vector<std::vector<double>> roir(src.getNrOfTimesteps());
		for(size_t j=0;j<src.getNrOfTimesteps();j++){
			data::TypedChunk<int> I=src.getChunk(0,0,s,j);
			if(j%2){ // odd
				Iodd+=I;
			} else {
				Ieven+=I;
			}
			Syt+= I*j; 
			Syy+= I*I;
			Stt+= j*j;
			
			S0++;
			St+=j;
			Stt+=j*j;

			//M:sub = I(X1:X2,Y1:Y2);
			data::MemChunk<int> sub(R,R);
			I.copyTileTo(sub,{X1,Y1});
			
			//M:roi(S0) = sum(sum(sub))/npx;
			roi.push_back(std::accumulate(sub.begin(),sub.end(),0.0)/(R*R));
			for(size_t r=0;r<R;r++){
				size_t x1 = xc - r/2;
				size_t y1 = yc - r/2; 

				data::MemChunk<int> sub(r+1,r+1);
				I.copyTileTo(sub,{x1,y1});
				roir[j].push_back(std::accumulate(sub.begin(),sub.end(),0.0)/(r*r));
			}
		}
		
		data::MemChunk<int> sub(R,R);
		(Iodd-Ieven).copyTileTo(sub,{X1,Y1});
		
		data::IOFactory::write(Iodd,"/tmp/odd.nii");
		data::IOFactory::write(Ieven,"/tmp/even.nii");
		data::IOFactory::write(Iodd-Ieven,"/tmp/diff.nii");

		//M:varI = var(sub(:));
		double mean=std::accumulate(sub.begin(),sub.end(),0.0)/(R*R);
		double varI=0;
		for(auto v:sub)
			varI+= (v-mean)*(v-mean)/(R*R);
		
		auto Iave = (Iodd+Ieven)* (1./src.getNrOfTimesteps());

	}

	
// 	ui->image_frame->layout()->addWidget(new isis::qt5::SimpleImageView(src));
// 
// 	auto bin=binarise(src);
// 	ui->bin_frame->layout()->addWidget(new isis::qt5::SimpleImageView(bin));
}

qualicheck::~qualicheck()
{
    delete ui;
}