1 + | |
2 + | |
3 + | |
4 + | |
5 + | |
6 + | |
7 + | |
8 + | |
9 + | |
10 + | |
11 + | |
12 + | |
13 + | |
14 + | |
15 + | |
16 + | |
17 + | |
18 + | |
19 + | |
20 + | |
21 + | |
22 + | |
23 + | |
24 + | |
25 + | |
26 + | |
27 + | |
28 + | |
29 + | |
30 + | |
31 + | |
32 + | |
33 + | |
34 + | |
35 + | #include <casacore/images/Images/ImageUtilities.h> |
36 + | #include <casa/iostream.h> |
37 + | #include <casa/aips.h> |
38 + | #include <casa/Exceptions/Error.h> |
39 + | #include <casa/BasicSL/String.h> |
40 + | #include <casa/Containers/Block.h> |
41 + | #include <measures/Measures/MRadialVelocity.h> |
42 + | #include <coordinates/Coordinates/CoordinateSystem.h> |
43 + | #include <casa/Logging/LogIO.h> |
44 + | #include <synthesis/ImagerObjects/SynthesisImager.h> |
45 + | #include <synthesis/ImagerObjects/SIImageStore.h> |
46 + | #include <imageanalysis/Utilities/SpectralImageUtil.h> |
47 + | #include <lattices/Lattices/LatticeConcat.h> |
48 + | #include <images/Images/PagedImage.h> |
49 + | #include <images/Images/ImageConcat.h> |
50 + | #include <images/Images/SubImage.h> |
51 + | #include <casa/namespace.h> |
52 + | #include <images/Images/TempImage.h> |
53 + | #include <coordinates/Coordinates/CoordinateUtil.h> |
54 + | #include <ms/MSSel/MSSourceIndex.h> |
55 + | #include <synthesis/ImagerObjects/SDMaskHandler.h> |
56 + | |
57 + | int main(int argc, char **argv) |
58 + | { |
59 + | using namespace std; |
60 + | using namespace casacore; |
61 + | using namespace casa; |
62 + | try{ |
63 + | |
64 + | |
65 + | String imagename; |
66 + | if (argc > 1) { |
67 + | imagename=argv[1]; |
68 + | } |
69 + | else { |
70 + | imagename="tLabelandFindRegions.temp.mask"; |
71 + | Int imsize = 100; |
72 + | Int nchan = 2; |
73 + | IPosition blc(4,45,50,0,0); |
74 + | IPosition trc(4,55,55,0,0); |
75 + | IPosition shape(4, imsize, imsize, 1, nchan); |
76 + | CoordinateSystem csys=CoordinateUtil::defaultCoords4D(); |
77 + | PagedImage<Float> maskImage(TiledShape(shape), csys, imagename); |
78 + | maskImage.setUnits(Unit("Jy/pixel")); |
79 + | maskImage.set(0.0); |
80 + | |
81 + | if (blc(0) <= trc(0) && blc(1) <= trc(1) && blc(3) < nchan) { |
82 + | Int dx = trc(0) - blc(0); |
83 + | Int dy = trc(1) - blc(1); |
84 + | for (uInt i=0+blc(3); i < trc(3)+1; i++) { |
85 + | for (uInt j=0; j < (uInt)dx+1; j++) { |
86 + | for (uInt k=0; k < (uInt)dy+1; k++) { |
87 + | IPosition loc(4,blc(0)+j, blc(1)+k, 0, Int(i)); |
88 + | Float val = 1.0; |
89 + | maskImage.putAt(val, loc); |
90 + | } |
91 + | } |
92 + | } |
93 + | } |
94 + | } |
95 + | |
96 + | CountedPtr<ImageInterface<Float> > im; |
97 + | im = ImageUtilities::openImage<Float>(imagename); |
98 + | IPosition shp=im->shape(); |
99 + | cerr<<"shp="<<shp<<endl; |
100 + | IPosition start(4, 0, 0, 0, 0); |
101 + | IPosition length(4, shp(0),shp(1),shp(2),1); |
102 + | Slicer sl(start, length); |
103 + | AxesSpecifier aspec(False); |
104 + | SubImage<Float>* subIm = new SubImage<Float>(*im, sl, aspec, True); |
105 + | cerr<<"subIm shape="<<subIm->shape()<<" slice="<<sl<<endl; |
106 + | SDMaskHandler smh; |
107 + | |
108 + | TempImage<Float>* tempIm = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates() ); |
109 + | tempIm->copyData(LatticeExpr<Float> (*subIm)); |
110 + | cerr<<"tempIm shape="<<tempIm->shape()<<endl; |
111 + | |
112 + | cerr<<"setup a blobmap ..."<<endl; |
113 + | TempImage<Float>* blobMap = new TempImage<Float> (TiledShape(subIm->shape(), subIm->niceCursorShape()), subIm->coordinates() ); |
114 + | blobMap->set(0); |
115 + | cerr<<"blobMap shape="<<blobMap->shape()<<endl; |
116 + | |
117 + | cerr<<"Calling labelRegions()..."<<endl; |
118 + | |
119 + | |
120 + | |
121 + | |
122 + | smh.labelRegions(*subIm, *blobMap); |
123 + | cerr<<"Calling findBlobSize()..."<<endl; |
124 + | Vector<Float> bloblist=smh.findBlobSize(*blobMap); |
125 + | for (uInt iblob=0; iblob<bloblist.nelements(); ++iblob) { |
126 + | cerr<<iblob<<" : "<<bloblist[iblob]<<endl; |
127 + | } |
128 + | |
129 + | }catch( AipsError e ){ |
130 + | cout << "Exception ocurred." << endl; |
131 + | cout << e.getMesg() << endl; |
132 + | } |
133 + | return 0; |
134 + | }; |