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