mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-04-19 18:40:01 +02:00
579 lines
15 KiB
C++
579 lines
15 KiB
C++
// SPDX-License-Identifier: LGPL-3.0-or-other
|
|
// Copyright (C) 2021 Contributors to the SLS Detector Package
|
|
#ifndef SLS_INTERPOLATION_H
|
|
#define SLS_INTERPOLATION_H
|
|
|
|
#include <cstdlib>
|
|
#ifndef MY_TIFF_IO_H
|
|
#include "tiffIO.h"
|
|
#endif
|
|
|
|
#ifndef DEF_QUAD
|
|
#define DEF_QUAD
|
|
enum quadrant {
|
|
TOP_LEFT=0,
|
|
TOP_RIGHT=1,
|
|
BOTTOM_LEFT=2,
|
|
BOTTOM_RIGHT=3,
|
|
UNDEFINED_QUADRANT=-1
|
|
};
|
|
#endif
|
|
#include <memory.h>
|
|
|
|
#include <stdio.h>
|
|
#include <iostream>
|
|
using namespace std;
|
|
|
|
//#ifdef MYROOT1
|
|
//: public TObject
|
|
//#endif
|
|
class slsInterpolation
|
|
{
|
|
|
|
public:
|
|
slsInterpolation(int nx=400, int ny=400, int ns=25, int nsy=-1) :nPixelsX(nx), nPixelsY(ny), nSubPixelsX(ns), nSubPixelsY(nsy),id(0) {
|
|
|
|
if (nSubPixelsY<=0) nSubPixelsY=nSubPixelsX;
|
|
hint=new int[nSubPixelsX*nx*nSubPixelsY*ny];
|
|
|
|
};
|
|
|
|
slsInterpolation(slsInterpolation *orig){
|
|
nPixelsX=orig->nPixelsX;
|
|
nPixelsY=orig->nPixelsY;
|
|
nSubPixelsX=orig->nSubPixelsX;
|
|
nSubPixelsY=orig->nSubPixelsY;
|
|
|
|
hint=new int[nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY];
|
|
memcpy(hint, orig->hint,nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY*sizeof(int));
|
|
|
|
};
|
|
|
|
virtual int setId(int i) {id=i; return id;};
|
|
|
|
virtual slsInterpolation* Clone() =0; /*{
|
|
return new slsInterpolation(this);
|
|
}*/
|
|
|
|
int getNSubPixelsX() {return nSubPixelsX;};
|
|
int getNSubPixelsY() {return nSubPixelsY;};
|
|
int getNSubPixels() {if (nSubPixelsX==nSubPixelsY) return nSubPixelsX; else return 0;};
|
|
void getNSubPixels(int &nsx, int &nsy) {nsx=nSubPixelsX; nsy=nsx=nSubPixelsY;}
|
|
|
|
void setNSubPixels(int ns, int nsy=-1) {
|
|
|
|
delete [] hint;
|
|
nSubPixelsX=ns;
|
|
if (nsy>0) nSubPixelsY=nsy;
|
|
else nSubPixelsY=ns;
|
|
|
|
hint=new int[nSubPixelsX*nPixelsX*nSubPixelsY*nPixelsY];
|
|
|
|
//return nSubPixels;
|
|
}
|
|
|
|
|
|
|
|
|
|
int getImageSize(int &nnx, int &nny, int &nsx, int &nsy) {
|
|
nnx=nSubPixelsX*nPixelsX;
|
|
nny=nSubPixelsY*nPixelsY;
|
|
nsx=nSubPixelsX;
|
|
nsy=nSubPixelsY;
|
|
return nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY;
|
|
};
|
|
|
|
|
|
int getImageSize() {
|
|
return nSubPixelsX*nSubPixelsY*nPixelsX*nPixelsY;
|
|
};
|
|
|
|
|
|
//create eta distribution, eta rebinnining etc.
|
|
//returns flat field image
|
|
virtual void prepareInterpolation(int &ok)=0;
|
|
|
|
//create interpolated image
|
|
//returns interpolated image
|
|
virtual int *getInterpolatedImage(){
|
|
// cout << "return interpolated image " << endl;
|
|
/* for (int i=0; i<nSubPixels* nSubPixels* nPixelsX*nPixelsY; i++) { */
|
|
/* cout << i << " " << hint[i] << endl; */
|
|
/* } */
|
|
return hint;
|
|
};
|
|
|
|
|
|
|
|
|
|
void *writeInterpolatedImage(const char * imgname) {
|
|
//cout << "!" <<endl;
|
|
float *gm=NULL;
|
|
int *dummy=getInterpolatedImage();
|
|
gm=new float[ nSubPixelsX* nSubPixelsY* nPixelsX*nPixelsY];
|
|
if (gm) {
|
|
for (int ix=0; ix<nPixelsX*nSubPixelsX; ix++) {
|
|
for (int iy=0; iy<nPixelsY*nSubPixelsY; iy++) {
|
|
gm[iy*nPixelsX*nSubPixelsX+ix]=dummy[iy*nPixelsX*nSubPixelsX+ix];
|
|
}
|
|
}
|
|
WriteToTiff(gm, imgname,nSubPixelsY* nPixelsX ,nSubPixelsY* nPixelsY);
|
|
delete [] gm;
|
|
} else cout << "Could not allocate float image " << endl;
|
|
return NULL;
|
|
}
|
|
|
|
|
|
|
|
//return position inside the pixel for the given photon
|
|
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0;
|
|
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0;
|
|
//return position inside the pixel for the given photon
|
|
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int quad, double &int_x, double &int_y)=0;
|
|
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cluster,double &etax, double &etay)=0;
|
|
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cluster,double &etax, double &etay)=0;
|
|
|
|
|
|
//return position inside the pixel for the given photon
|
|
virtual void clearInterpolatedImage() {
|
|
|
|
|
|
for (int ix=0; ix<nPixelsX*nSubPixelsX; ix++) {
|
|
for (int iy=0; iy<nPixelsY*nSubPixelsY; iy++) {
|
|
hint[iy*nPixelsX*nSubPixelsX+ix]=0;
|
|
}
|
|
}
|
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
virtual int *addToImage(double int_x, double int_y){
|
|
int iy=((double)nSubPixelsY)*int_y;
|
|
int ix=((double)nSubPixelsX)*int_x;
|
|
if (ix>=0 && ix<(nPixelsX*nSubPixelsX) && iy<(nSubPixelsY*nPixelsY) && iy>=0 ){
|
|
// cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << " " << hint[ix+iy*nPixelsX*nSubPixels];
|
|
(*(hint+ix+iy*nPixelsX*nSubPixelsX))+=1;
|
|
// cout << " " << hint[ix+iy*nPixelsX*nSubPixels] << endl;
|
|
}// else
|
|
// cout << "bad! "<< int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl;
|
|
|
|
return hint;
|
|
};
|
|
|
|
|
|
virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
|
|
virtual int addToFlatField(int *cluster, double &etax, double &etay)=0;
|
|
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0;
|
|
virtual int addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0;
|
|
virtual int addToFlatField(double etax, double etay)=0;
|
|
|
|
virtual int *getFlatField(){return NULL;};
|
|
virtual int *setFlatField(int *h, int nb=-1, double emin=-1, double emax=-1){return NULL;};
|
|
virtual void *writeFlatField(const char * imgname){return NULL;};
|
|
virtual void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){return NULL;};
|
|
virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();};
|
|
|
|
virtual void resetFlatField()=0;
|
|
|
|
//virtual void Streamer(TBuffer &b);
|
|
|
|
static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){
|
|
double cli[3*3];//=new int[3*3];
|
|
for (int i=0; i<9; i++)
|
|
cli[i]=cl[i];
|
|
return calcQuad(cli, sum, totquad, sDum);
|
|
|
|
}
|
|
|
|
|
|
static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
|
|
|
|
int corner = UNDEFINED_QUADRANT;
|
|
/* double *cluster[3]; */
|
|
/* cluster[0]=cl; */
|
|
/* cluster[1]=cl+3; */
|
|
/* cluster[2]=cl+6; */
|
|
|
|
sum=0;
|
|
int xoff=0, yoff=0;
|
|
#ifndef WRITE_QUAD
|
|
double sumBL=0;
|
|
double sumTL=0;
|
|
double sumBR=0;
|
|
double sumTR=0;
|
|
for (int ix=0; ix<3; ix++) {
|
|
for (int iy=0; iy<3; iy++) {
|
|
sum+=cl[ix+3*iy];
|
|
if (ix<=1 && iy<=1) sumBL+=cl[ix+iy*3];
|
|
if (ix<=1 && iy>=1) sumTL+=cl[ix+iy*3];
|
|
if (ix>=1 && iy<=1) sumBR+=cl[ix+iy*3];
|
|
if (ix>=1 && iy>=1) sumTR+=cl[ix+iy*3];
|
|
}
|
|
}
|
|
|
|
/* sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0]; */
|
|
/* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */
|
|
corner = BOTTOM_LEFT;
|
|
totquad=sumBL;
|
|
xoff=0;
|
|
yoff=0;
|
|
|
|
|
|
if(sumTL >= totquad){
|
|
/* #ifdef WRITE_QUAD */
|
|
/* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
|
|
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
|
|
/* if (sumTL ==sum) { */
|
|
/* #endif */
|
|
corner = TOP_LEFT;
|
|
totquad=sumTL;
|
|
xoff=0;
|
|
yoff=1;
|
|
/* #ifdef WRITE_QUAD */
|
|
/* } */
|
|
/* #endif */
|
|
}
|
|
|
|
if(sumBR >= totquad){
|
|
/* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */
|
|
/* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */
|
|
|
|
/* #ifdef WRITE_QUAD */
|
|
/* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
|
|
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
|
|
/* if (sumBR ==sum) { */
|
|
/* #endif */
|
|
corner = BOTTOM_RIGHT;
|
|
xoff=1;
|
|
yoff=0;
|
|
totquad=sumBR;
|
|
/* #ifdef WRITE_QUAD */
|
|
/* } */
|
|
/* #endif */
|
|
}
|
|
|
|
if(sumTR >= totquad){
|
|
/* #ifdef WRITE_QUAD */
|
|
/* /\* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; *\/ */
|
|
/* /\* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; *\/ */
|
|
/* if (sumTR ==sum) { */
|
|
/* #endif */
|
|
xoff=1;
|
|
yoff=1;
|
|
/* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */
|
|
/* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */
|
|
corner = TOP_RIGHT;
|
|
totquad=sumTR;
|
|
/* #ifdef WRITE_QUAD */
|
|
/* } */
|
|
/* #endif */
|
|
}
|
|
|
|
#endif
|
|
#ifdef WRITE_QUAD
|
|
|
|
double sumB=0;
|
|
double sumT=0;
|
|
double sumR=0;
|
|
double sumL=0;
|
|
|
|
for (int ix=0; ix<3; ix++) {
|
|
for (int iy=0; iy<3; iy++) {
|
|
sum+=cl[ix+3*iy];
|
|
if (ix<1 ) sumL+=cl[ix+iy*3];
|
|
if (ix>1) sumR+=cl[ix+iy*3];
|
|
if (iy<1) sumB=cl[ix+iy*3];
|
|
if (iy>1) sumT+=cl[ix+iy*3];
|
|
}
|
|
}
|
|
|
|
totquad=sum;
|
|
if ( sumT==0 && sumR==0) {
|
|
corner = BOTTOM_LEFT;
|
|
xoff=0;
|
|
yoff=0;
|
|
} else if ( sumB==0 && sumR==0 ) {
|
|
corner = TOP_LEFT;
|
|
xoff=0;
|
|
yoff=1;
|
|
} else if ( sumT==0 && sumL==0) {
|
|
corner = BOTTOM_RIGHT;
|
|
xoff=1;
|
|
yoff=0;
|
|
} else if ( sumB==0 && sumL==0) {
|
|
xoff=1;
|
|
yoff=1;
|
|
corner = TOP_RIGHT;
|
|
} else
|
|
printf("** bad 2x2 cluster!\n");
|
|
|
|
#endif
|
|
|
|
|
|
for (int ix=0; ix<2; ix++) {
|
|
for (int iy=0; iy<2; iy++) {
|
|
sDum[iy][ix] = cl[ix+xoff+(iy+yoff)*3];
|
|
}
|
|
}
|
|
|
|
return corner;
|
|
|
|
}
|
|
|
|
static int calcEta(double totquad, double sDum[2][2], double &etax, double &etay){
|
|
double t,r;
|
|
|
|
if (totquad>0) {
|
|
t = sDum[1][0] + sDum[1][1];
|
|
r = sDum[0][1] + sDum[1][1];
|
|
etax=r/totquad;
|
|
etay=t/totquad;
|
|
}
|
|
return 0;
|
|
|
|
}
|
|
|
|
static int calcEta(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEta(totquad, sDum, etax, etay);
|
|
|
|
return corner;
|
|
}
|
|
|
|
|
|
static int calcEta(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEta(totquad, sDum, etax, etay);
|
|
|
|
return corner;
|
|
}
|
|
|
|
|
|
static int calcEtaL(double totquad, int corner, double sDum[2][2], double &etax, double &etay){
|
|
double t,r, toth, totv;
|
|
if (totquad>0) {
|
|
switch(corner) {
|
|
case TOP_LEFT:
|
|
t = sDum[1][1];
|
|
r = sDum[0][1] ;
|
|
toth=sDum[0][1]+sDum[0][0];
|
|
totv=sDum[0][1]+sDum[1][1];
|
|
break;
|
|
case TOP_RIGHT:
|
|
t = sDum[1][0] ;
|
|
r = sDum[0][1] ;
|
|
toth=sDum[0][1]+sDum[0][0];
|
|
totv=sDum[1][0]+sDum[0][0];
|
|
break;
|
|
case BOTTOM_LEFT:
|
|
r = sDum[1][1] ;
|
|
t = sDum[1][1] ;
|
|
toth=sDum[1][0]+sDum[1][1];
|
|
totv=sDum[0][1]+sDum[1][1];
|
|
break;
|
|
case BOTTOM_RIGHT:
|
|
t = sDum[1][0] ;
|
|
r = sDum[1][1] ;
|
|
toth=sDum[1][0]+sDum[1][1];
|
|
totv=sDum[1][0]+sDum[0][0];
|
|
break;
|
|
default:
|
|
etax=-1000;
|
|
etay=-1000;
|
|
return 0;
|
|
}
|
|
//etax=r/totquad;
|
|
//etay=t/totquad;
|
|
etax=r/toth;
|
|
etay=t/totv;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
static int calcEtaL(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEtaL(totquad, corner, sDum, etax, etay);
|
|
|
|
return corner;
|
|
}
|
|
|
|
static int calcEtaL(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEtaL(totquad, corner, sDum, etax, etay);
|
|
|
|
return corner;
|
|
}
|
|
|
|
|
|
|
|
static int calcEtaC3(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
|
|
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEta(sum, sDum, etax, etay);
|
|
return corner;
|
|
|
|
}
|
|
|
|
|
|
|
|
static int calcEtaC3(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
|
|
|
|
int corner = calcQuad(cl,sum,totquad,sDum);
|
|
calcEta(sum, sDum, etax, etay);
|
|
return corner;
|
|
|
|
}
|
|
|
|
|
|
|
|
static int calcEta3(double *cl, double &etax, double &etay, double &sum) {
|
|
double l=0,r=0,t=0,b=0, val;
|
|
sum=0;
|
|
// int quad;
|
|
for (int ix=0; ix<3; ix++) {
|
|
for (int iy=0; iy<3; iy++) {
|
|
val=cl[ix+3*iy];
|
|
sum+=val;
|
|
if (iy==0) l+=val;
|
|
if (iy==2) r+=val;
|
|
if (ix==0) b+=val;
|
|
if (ix==2) t+=val;
|
|
}
|
|
}
|
|
if (sum>0) {
|
|
etax=(-l+r)/sum;
|
|
etay=(-b+t)/sum;
|
|
}
|
|
|
|
if (etax>=0 && etay>=0)
|
|
return TOP_RIGHT;
|
|
if (etax<0 && etay>=0)
|
|
return TOP_LEFT;
|
|
if (etax<0 && etay<0)
|
|
return BOTTOM_LEFT;
|
|
return BOTTOM_RIGHT;
|
|
}
|
|
|
|
|
|
static int calcEta3(int *cl, double &etax, double &etay, double &sum) {
|
|
double cli[9];
|
|
for (int ix=0; ix<9; ix++) cli[ix]=cl[ix];
|
|
|
|
return calcEta3(cli, etax, etay, sum);
|
|
}
|
|
|
|
|
|
|
|
static int calcEta3X(double *cl, double &etax, double &etay, double &sum) {
|
|
double l,r,t,b;
|
|
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
|
|
if (sum>0) {
|
|
l=cl[3];
|
|
r=cl[5];
|
|
b=cl[1];
|
|
t=cl[7];
|
|
etax=(-l+r)/sum;
|
|
etay=(-b+t)/sum;
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
static int calcEta3X(int *cl, double &etax, double &etay, double &sum) {
|
|
double cli[9];
|
|
for (int ix=0; ix<9; ix++) cli[ix]=cl[ix];
|
|
|
|
return calcEta3X(cli, etax, etay, sum);
|
|
}
|
|
|
|
|
|
/* static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { */
|
|
/* double l,r,t,b, sum; */
|
|
/* int yoff; */
|
|
/* switch (quad) { */
|
|
/* case BOTTOM_LEFT: */
|
|
/* case BOTTOM_RIGHT: */
|
|
/* yoff=0; */
|
|
/* break; */
|
|
/* case TOP_LEFT: */
|
|
/* case TOP_RIGHT: */
|
|
/* yoff=1; */
|
|
/* break; */
|
|
/* default: */
|
|
/* ; */
|
|
/* } */
|
|
/* l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
|
|
/* r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
|
|
/* b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
|
|
/* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
|
|
/* sum=t+b; */
|
|
/* if (sum>0) { */
|
|
/* etax=(-l+r)/sum; */
|
|
/* etay=(+t)/sum; */
|
|
/* } */
|
|
|
|
/* return -1; */
|
|
/* } */
|
|
|
|
/* static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { */
|
|
/* double l,r,t,b, sum; */
|
|
/* int yoff; */
|
|
/* switch (quad) { */
|
|
/* case BOTTOM_LEFT: */
|
|
/* case BOTTOM_RIGHT: */
|
|
/* yoff=0; */
|
|
/* break; */
|
|
/* case TOP_LEFT: */
|
|
/* case TOP_RIGHT: */
|
|
/* yoff=1; */
|
|
/* break; */
|
|
/* default: */
|
|
/* ; */
|
|
/* } */
|
|
/* l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
|
|
/* r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
|
|
/* b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
|
|
/* t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
|
|
/* sum=t+b; */
|
|
/* if (sum>0) { */
|
|
/* etax=(-l+r)/sum; */
|
|
/* etay=(+t)/sum; */
|
|
/* } */
|
|
|
|
/* return -1; */
|
|
/* } */
|
|
|
|
|
|
|
|
/* static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { */
|
|
/* double l,r,t,b; */
|
|
/* sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8]; */
|
|
/* if (sum>0) { */
|
|
/* l=cl[3]; */
|
|
/* r=cl[5]; */
|
|
/* b=cl[1]; */
|
|
/* t=cl[7]; */
|
|
/* etax=(-l+r)/sum; */
|
|
/* etay=(-b+t)/sum; */
|
|
/* } */
|
|
/* return -1; */
|
|
/* } */
|
|
|
|
|
|
|
|
|
|
|
|
|
|
protected:
|
|
int nPixelsX, nPixelsY;
|
|
int nSubPixelsX, nSubPixelsY;
|
|
int id;
|
|
int *hint;
|
|
};
|
|
|
|
#endif
|