Interpolation in both directions for rect pixels works, saving of the rectangular image is done by rebinning

This commit is contained in:
2020-01-15 17:27:46 +01:00
parent 27f2321d64
commit e4ab4547db
4 changed files with 196 additions and 86 deletions

View File

@ -196,11 +196,12 @@ class slsInterpolation
/* 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;
int xoff=0, yoff=0;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
sum+=cl[ix+3*iy];
@ -220,34 +221,95 @@ class slsInterpolation
if(sumTL >= totquad){
/* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; */
/* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; */
/* #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;
corner = BOTTOM_RIGHT;
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++) {

View File

@ -12,7 +12,7 @@ MAIN=moench03ClusterFinder.cpp
#-lhdf5
#DESTDIR?=../bin
all: moenchClusterFinderRect moenchMakeEtaRect moenchInterpolationRect moenchNoInterpolationRect moenchPhotonCounterRect moenchAnalogRect
all: moenchClusterFinderRect moenchMakeEtaRect moenchInterpolationRect moenchNoInterpolationRect moenchPhotonCounterRect moenchAnalogRect moenchMakeEtaRect1 moenchInterpolationRect1
@ -34,6 +34,12 @@ moenchPhotonCounterRect: moenchPhotonCounter.cpp $(INCS) clean
moenchAnalogRect: moenchPhotonCounter.cpp $(INCS) clean
g++ -o moenchAnalogRect moenchPhotonCounter.cpp $(LDFLAG) $(INCDIR) -DNEWRECEIVER -DANALOG -DRECT
moenchMakeEtaRect1: moench03Interpolation.cpp $(INCS) clean
g++ -o moenchMakeEtaRect1 moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DFF -DRECT -DRECT1
moenchInterpolationRect1: moench03Interpolation.cpp $(INCS) clean
g++ -o moenchInterpolationRect1 moench03Interpolation.cpp $(LDFLAG) $(INCDIR) -DRECT -DRECT1
clean:
rm -f moenchClusterFinderRect moenchMakeEtaRect moenchInterpolationRect moenchNoInterpolationRect moenchPhotonCounterRect moenchAnalogRect

View File

@ -5,7 +5,7 @@
//#include "moench03T1ZmqData.h"
//#define DOUBLE_SPH
//#define MANYFILES
#define WRITE_QUAD
#ifdef DOUBLE_SPH
#include "single_photon_hit_double.h"
#endif
@ -101,10 +101,14 @@ int main(int argc, char *argv[]) {
#endif
int nSubPixels=nsubpix;
int nSubPixelsY=nsubpix;
#ifdef RECT
nSubPixels=1; //might make more sense using 2?
etabins=2;
#ifdef RECT1
// nSubPixels=2; //might make more sense using 2?
// etabins=100;
float *rectimg=new float[NC*NR*nSubPixelsY];;
#endif
// #endif
#ifndef NOINTERPOLATION
eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nSubPixels, nSubPixelsY, etabins, etabinsY, etamin, etamax);
//eta2InterpolationCleverAdaptiveBins *interp=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax);
@ -172,6 +176,8 @@ int main(int argc, char *argv[]) {
// #ifdef SOLEIL
// if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) {
// #endif
if (cl.x<0 || cl.y<0 || cl.x>NC || cl.y>NR) {
} else {
#ifndef FF
// interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);
interp->getInterpolatedPosition(cl.x,cl.y, etax, etay, quad,int_x, int_y);
@ -181,14 +187,14 @@ int main(int argc, char *argv[]) {
// cout << int_x << " " << int_y << endl;
// cout <<"**************"<< endl;
// if (etax!=0 && etay!=0 && etax!=1 && etay!=1)
interp->addToImage(int_x, int_y);
if (int_x<0 || int_y<0 || int_x>NC || int_y>NR) {
if (int_x<-1 || int_y<-1 || int_x>NC+1 || int_y>NR+1) {
cout <<"**************"<< endl;
cout << cl.x << " " << cl.y << " " << sum << endl;
cl.print();
cout << int_x << " " << int_y << endl;
cout <<"**************"<< endl;
}
} else
interp->addToImage(int_x, int_y);
#endif
#ifdef FF
// interp->addToFlatField(cl.get_cluster(), etax, etay);
@ -200,6 +206,7 @@ int main(int argc, char *argv[]) {
// if (etax==0 || etay==0) cout << cl.x << " " << cl.y << endl;
#endif
}
// #ifdef SOLEIL
// }
// #endif
@ -225,18 +232,31 @@ int main(int argc, char *argv[]) {
#endif
#ifndef FF
#ifndef RECT1
interp->writeInterpolatedImage(outfname);
#endif
img=interp->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nSubPixels; isx++) {
for (isy=0; isy<nSubPixelsY; isy++) {
totimg[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)]+=img[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)];
#ifdef RECT1
if (isx==0)
rectimg[ix+(iy*nSubPixelsY+isy)*(NC)]=img[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)];
else
rectimg[ix+(iy*nSubPixelsY+isy)*(NC)]+=img[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)];
#endif
}
}
}
}
#ifdef RECT1
WriteToTiff(rectimg, outfname,NC,NR*nSubPixelsY);
#endif
cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ") nph="<< nph <<endl;
interp->clearInterpolatedImage();
#endif
@ -246,8 +266,30 @@ int main(int argc, char *argv[]) {
}
#ifndef FF
sprintf(outfname,argv[3],11111);
#ifndef RECT1
WriteToTiff(totimg, outfname,NC*nSubPixels,NR*nSubPixelsY);
#endif
#ifdef RECT1
img=interp->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nSubPixels; isx++) {
for (isy=0; isy<nSubPixelsY; isy++) {
if (isx==0)
rectimg[ix+(iy*nSubPixelsY+isy)*(NC)]=totimg[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)];
else
rectimg[ix+(iy*nSubPixelsY+isy)*(NC)]+=totimg[ix*nSubPixels+isx+(iy*nSubPixelsY+isy)*(NC*nSubPixels)];
}
}
}
}
WriteToTiff(rectimg, outfname,NC,NR*nSubPixelsY);
#endif
#endif
#ifdef FF
interp->writeFlatField(outfname);

View File

@ -119,44 +119,44 @@ class single_photon_hit {
if (fread((void*)qq, 1, 4*sizeof(int), myFile)) {
quad=TOP_RIGHT;
/* int mm=qq[0]; */
/* for (int i=1; i<4; i++) { */
/* if (qq[i]>mm) { */
/* switch (i) { */
/* case 1: */
/* quad=TOP_LEFT; */
/* break; */
/* case 2: */
/* quad=BOTTOM_RIGHT; */
/* break; */
/* case 3: */
/* quad=BOTTOM_LEFT; */
/* break; */
/* default: */
/* ; */
/* } */
int mm=qq[0];
for (int i=1; i<4; i++) {
if (qq[i]>mm) {
switch (i) {
case 1:
quad=TOP_LEFT;
break;
case 2:
quad=BOTTOM_RIGHT;
break;
case 3:
quad=BOTTOM_LEFT;
break;
default:
;
}
mm=qq[i];
}
/* } */
/* } */
}
/* switch(quad) { */
/* case TOP_LEFT: */
/* data[0]=0; */
/* data[1]=0; */
/* data[2]=0; */
/* data[3]=qq[0]; */
/* data[4]=qq[1]; */
/* data[5]=0; */
/* data[6]=qq[2]; */
/* data[7]=qq[3]; */
/* data[8]=0; */
/* x=x+1; */
/* y=y; */
/* break; */
switch(quad) {
case TOP_LEFT:
data[0]=0;
data[1]=0;
data[2]=0;
data[3]=qq[0];
data[4]=qq[1];
data[5]=0;
data[6]=qq[2];
data[7]=qq[3];
data[8]=0;
x=x+1;
y=y;
break;
/* case TOP_RIGHT: */
case TOP_RIGHT:
data[0]=0;
data[1]=0;
data[2]=0;
@ -168,40 +168,40 @@ class single_photon_hit {
data[8]=qq[3];
x=x;
y=y;
/* break; */
break;
/* case BOTTOM_LEFT: */
/* data[0]=qq[0]; */
/* data[1]=qq[1]; */
/* data[2]=0; */
/* data[3]=qq[2]; */
/* data[4]=qq[3]; */
/* data[5]=0; */
/* data[6]=0; */
/* data[7]=0; */
/* data[8]=0; */
/* x=x+1; */
/* y=y+1; */
/* break; */
/* case BOTTOM_RIGHT: */
/* data[0]=0; */
/* data[1]=qq[0]; */
/* data[2]=qq[1]; */
/* data[3]=0; */
/* data[4]=qq[2]; */
/* data[5]=qq[3]; */
/* data[6]=0; */
/* data[7]=0; */
/* data[8]=0; */
/* x=x; */
/* y=y+1; */
/* break; */
case BOTTOM_LEFT:
data[0]=qq[0];
data[1]=qq[1];
data[2]=0;
data[3]=qq[2];
data[4]=qq[3];
data[5]=0;
data[6]=0;
data[7]=0;
data[8]=0;
x=x+1;
y=y+1;
break;
case BOTTOM_RIGHT:
data[0]=0;
data[1]=qq[0];
data[2]=qq[1];
data[3]=0;
data[4]=qq[2];
data[5]=qq[3];
data[6]=0;
data[7]=0;
data[8]=0;
x=x;
y=y+1;
break;
/* default: */
/* ; */
/* } */
default:
;
}
return 1;
}
#endif