1 | #include <iostream> |
---|
2 | #include <iomanip> |
---|
3 | //#include <Cubes/cubes.hh> |
---|
4 | #include <Utils/utils.hh> |
---|
5 | |
---|
6 | using namespace std; |
---|
7 | |
---|
8 | void printDetails(int num,Detection &obj); |
---|
9 | DetectionList *search3DBits(Cube &cube, Param &par) |
---|
10 | { |
---|
11 | DetectionList *objectList = new DetectionList; |
---|
12 | int count=0; |
---|
13 | float nullPixValue = par.getNullPixVal(); |
---|
14 | bool flagFDR = par.getFlagFDR(); |
---|
15 | |
---|
16 | for(int bit=0; bit<3; bit++){ |
---|
17 | |
---|
18 | long xdim = cube.getDimX(); |
---|
19 | long ydim = cube.getDimY(); |
---|
20 | long zdim = cube.getDimZ() / 2; |
---|
21 | long size = xdim*ydim*zdim; |
---|
22 | |
---|
23 | float *arr = new float[size]; |
---|
24 | float *newarr = new float[size]; |
---|
25 | int goodSize=0; |
---|
26 | for(int x=0;x<xdim;x++){ |
---|
27 | for(int y=0;y<ydim;y++){ |
---|
28 | for(int z=0;z<zdim;z++){ |
---|
29 | long arrpos = x + y*xdim + z*xdim*ydim; |
---|
30 | long cubepos = x + y*xdim + (z + bit*zdim/2)*xdim*ydim; |
---|
31 | arr[arrpos] = cube.getPixValue(cubepos); |
---|
32 | if((!par.getFlagNullPix())||(cube.getPixValue(cubepos)!=nullPixValue)) |
---|
33 | newarr[goodSize++] = cube.getPixValue(cubepos); |
---|
34 | } |
---|
35 | } |
---|
36 | } |
---|
37 | |
---|
38 | float mean,sigma; |
---|
39 | if(goodSize>0){ |
---|
40 | findMedianStats(newarr,goodSize,mean,sigma); |
---|
41 | cerr<<xdim<<"x"<<ydim<<"x"<<zdim<<"-->"; |
---|
42 | cerr<<par.getCut()<<" "<<par.getAtrousCut()<<" "<<par.getFlagFDR()<<endl; |
---|
43 | atrous3DReconstruct(xdim,ydim,zdim,arr,newarr,par); |
---|
44 | //atrous3DReconstructBetter(xdim,ydim,zdim,arr,newarr,par); |
---|
45 | cerr<<"<--"<<endl; |
---|
46 | |
---|
47 | // cube.saveArray(newarr,size); |
---|
48 | // cube.saveOrig(arr,size); |
---|
49 | |
---|
50 | /* Search spatially */ |
---|
51 | long dim[2]={xdim,ydim}; |
---|
52 | for(int channel=0;channel<zdim;channel++){ |
---|
53 | bool skipThisOne = par.getFlagMW() && |
---|
54 | (channel>=par.getMinMW()) && (channel<=par.getMaxMW()); // ignore the Milky Way. |
---|
55 | if(!skipThisOne){ |
---|
56 | Image *slice = new Image(dim); |
---|
57 | for(long j=0;j<xdim*ydim;j++) |
---|
58 | slice->setPixValue( j, newarr[j+channel*xdim*ydim]); |
---|
59 | |
---|
60 | float *temparr = new float[xdim*ydim]; |
---|
61 | int ct=0; |
---|
62 | for(int i=0;i<xdim*ydim;i++) |
---|
63 | if((!par.getFlagNullPix())||(arr[channel*xdim*ydim+i]!=nullPixValue)) |
---|
64 | temparr[ct++] = arr[channel*xdim*ydim+i]; |
---|
65 | findMedianStats(temparr,ct,mean,sigma); |
---|
66 | delete [] temparr; |
---|
67 | |
---|
68 | slice->setStats(mean,sigma,par.getCut()); |
---|
69 | if(flagFDR) setupFDR(*slice,par); |
---|
70 | lutz_detect(*slice,par); |
---|
71 | for(int i=0;i<slice->getNumObj();i++){ |
---|
72 | count++; |
---|
73 | Detection *obj = new Detection; |
---|
74 | *obj = slice->getObject(i); |
---|
75 | for(int j=0;j<obj->getSize();j++) obj->setZ(j,channel); |
---|
76 | obj->calcParams(); |
---|
77 | objectList->addObject(*obj); |
---|
78 | printDetails(count,*obj); |
---|
79 | delete obj; |
---|
80 | } |
---|
81 | for(int i=0;i<slice->getSize();i++) |
---|
82 | cube.setDetectMapValue(i,cube.getDetectMapValue(i) + slice->getMaskValue(i)); |
---|
83 | delete slice; |
---|
84 | } |
---|
85 | } |
---|
86 | |
---|
87 | /* Search in velocity */ |
---|
88 | cerr<<endl; |
---|
89 | long dim2[2]={zdim,1}; |
---|
90 | for(int npixel=0;npixel<xdim*ydim;npixel++){ |
---|
91 | Image *spectrum = new Image(dim2); |
---|
92 | for(long zpos=0;zpos<zdim;zpos++) |
---|
93 | spectrum->setPixValue( zpos, newarr[npixel+zpos*xdim*ydim] ); |
---|
94 | |
---|
95 | float *temparr = new float[zdim]; |
---|
96 | int ct=0; |
---|
97 | for(int i=0;i<zdim;i++) |
---|
98 | if((!par.getFlagNullPix())||(arr[npixel+i*xdim*ydim]!=nullPixValue)) |
---|
99 | temparr[ct++] = arr[npixel+i*xdim*ydim]; |
---|
100 | findMedianStats(temparr,zdim,mean,sigma); |
---|
101 | delete [] temparr; |
---|
102 | |
---|
103 | spectrum->setStats(mean,sigma,par.getCut()); |
---|
104 | if(flagFDR) setupFDR(*spectrum,par); |
---|
105 | lutz_detect(*spectrum,par); |
---|
106 | for(int i=0;i<spectrum->getNumObj();i++){ |
---|
107 | count++; |
---|
108 | Detection *obj = new Detection; |
---|
109 | *obj = spectrum->getObject(i); |
---|
110 | for(int j=0;j<obj->getSize();j++){ |
---|
111 | obj->setZ(j,obj->getX(j)); |
---|
112 | obj->setX(j,npixel%xdim); |
---|
113 | obj->setY(j,npixel/xdim); |
---|
114 | } |
---|
115 | obj->calcParams(); |
---|
116 | objectList->addObject(*obj); |
---|
117 | printDetails(count,*obj); |
---|
118 | delete obj; |
---|
119 | } |
---|
120 | for(int i=0;i<spectrum->getSize();i++) |
---|
121 | cube.setDetectMapValue(npixel,cube.getDetectMapValue(npixel)+spectrum->getMaskValue(i)); |
---|
122 | delete spectrum; |
---|
123 | } |
---|
124 | |
---|
125 | } |
---|
126 | delete [] arr; |
---|
127 | delete [] newarr; |
---|
128 | } |
---|
129 | |
---|
130 | cerr<<"Found "<<count<<" objects.\n"; |
---|
131 | return objectList; |
---|
132 | } |
---|
133 | |
---|
134 | // void printDetails(int num,Detection &obj) |
---|
135 | // { |
---|
136 | // cerr<<"Obj#"; |
---|
137 | // cerr<<setw(4)<<num<<": "; |
---|
138 | // cerr.setf(ios::fixed); |
---|
139 | // cerr<<setw(5)<<setprecision(1)<<obj.getXcentre()<<" "; |
---|
140 | // cerr<<setw(5)<<setprecision(1)<<obj.getYcentre()<<" "; |
---|
141 | // cerr<<setw(3)<<setprecision(1)<<obj.getZcentre()<<" "; |
---|
142 | // cerr<<" npix = "<<obj.getSize(); |
---|
143 | // cerr<<" F = "<<obj.getTotalFlux()<<endl; |
---|
144 | // } |
---|