Skip to content

Commit 425fc8b

Browse files
bourgeslbourgesl
committed
process both OI_VIS and OI_VIS2 tables but filter out redudant spatial coordinates after very small rounding (~ precision) on (U,V, lambda) tuples
1 parent 97023d5 commit 425fc8b

File tree

1 file changed

+180
-89
lines changed

1 file changed

+180
-89
lines changed

‎src/main/java/fr/jmmc/oitools/processing/BeamEstimator.java

+180-89
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import fr.jmmc.oitools.model.OIData;
88
import fr.jmmc.oitools.model.OIFitsFile;
99
import fr.jmmc.oitools.model.OIFitsLoader;
10+
import fr.jmmc.oitools.model.OIVis;
1011
import fr.jmmc.oitools.model.OIVis2;
1112
import java.awt.BasicStroke;
1213
import java.awt.Color;
@@ -17,6 +18,7 @@
1718
import java.util.Arrays;
1819
import java.util.Collection;
1920
import java.util.Locale;
21+
import java.util.TreeSet;
2022
import java.util.logging.Level;
2123
import javax.swing.JFrame;
2224
import javax.swing.JPanel;
@@ -40,107 +42,123 @@ public class BeamEstimator {
4042
private final static double SCALE_FACTOR = STDDEV_TO_HWHM / (2.0 * Math.PI);
4143
/** Specify the value of one milli arcsecond in degrees */
4244
public static final double DEG_IN_MILLI_ARCSEC = 3600000d;
45+
/** rounding precision on (U,V) coordinates in meter ~ 5 mm */
46+
public static final double PREC_UV = 5e-3;
47+
/** rounding precision on wavelengths in meter ~ 5e-11 m = 0.5 angstrom */
48+
public static final double PREC_WL = 5e-11;
4349

4450
public static BeamInfo computeBeamInfo(final SelectorResult result) {
45-
int n = 0;
51+
int n = 0, t = 0;
4652
double s11 = 0.0;
4753
double s22 = 0.0;
4854
double s12 = 0.0;
4955

50-
// TODO: how to handle only OIVIS or OIT3 coords or mix (VIS/VIS2/T3) tables that may have repeated U-V coords ?
56+
final TreeSet<UVTuple> uvSet = new TreeSet<>();
57+
5158
for (final OIData oiData : result.getOIDatas()) {
52-
final double[][] u; // rad-1
53-
final double[][] v; // rad-1
59+
final double[] ucoord; // m
60+
final double[] vcoord; // m
5461

5562
if (oiData instanceof OIVis2) {
5663
final OIVis2 vis2 = (OIVis2) oiData;
57-
58-
u = vis2.getSpatialUCoord();
59-
v = vis2.getSpatialVCoord();
60-
/*
61-
}
62-
else if (oiData instanceof OIVis) {
64+
ucoord = vis2.getUCoord();
65+
vcoord = vis2.getVCoord();
66+
} else if (oiData instanceof OIVis) {
6367
final OIVis vis = (OIVis) oiData;
64-
65-
u = vis.getSpatialUCoord();
66-
v = vis.getSpatialVCoord();
67-
*/
68+
ucoord = vis.getUCoord();
69+
vcoord = vis.getVCoord();
6870
} else {
71+
/** ignore OI_T3 as their UV coordinates should be redudant with OI_VIS/OI_VIS2 tables by design */
6972
continue;
7073
}
7174

72-
// get the optional masks for this OIData table:
73-
final IndexMask maskOIData1D = result.getDataMask1DNotFull(oiData);
74-
final IndexMask maskOIData2D = result.getDataMask2DNotFull(oiData);
75-
// get the optional wavelength mask for the OIData's wavelength table:
76-
final IndexMask maskWavelength = result.getWavelengthMaskNotFull(oiData.getOiWavelength());
77-
78-
final int idxNone = (maskOIData2D != null) ? maskOIData2D.getIndexNone() : -1;
79-
final int idxFull = (maskOIData2D != null) ? maskOIData2D.getIndexFull() : -1;
80-
8175
final int nRows = oiData.getNbRows();
8276
final int nWaves = oiData.getNWave();
8377

84-
final boolean[][] flags = oiData.getFlag();
85-
86-
IndexMask maskOIData2DRow = null;
78+
if (nWaves != 0) {
79+
final double[] effWaves = oiData.getOiWavelength().getEffWaveAsDouble();
8780

88-
// Iterate on table rows (i):
89-
for (int i = 0; i < nRows; i++) {
81+
// get the optional masks for this OIData table:
82+
final IndexMask maskOIData1D = result.getDataMask1DNotFull(oiData);
83+
final IndexMask maskOIData2D = result.getDataMask2DNotFull(oiData);
84+
// get the optional wavelength mask for the OIData's wavelength table:
85+
final IndexMask maskWavelength = result.getWavelengthMaskNotFull(oiData.getOiWavelength());
9086

91-
// check optional data mask 1D:
92-
if ((maskOIData1D != null) && !maskOIData1D.accept(i)) {
93-
// if bit is false for this row, we hide this row
94-
continue;
95-
}
87+
final int idxNone = (maskOIData2D != null) ? maskOIData2D.getIndexNone() : -1;
88+
final int idxFull = (maskOIData2D != null) ? maskOIData2D.getIndexFull() : -1;
9689

97-
// check mask 2D for row None flag:
98-
if (maskOIData2D != null) {
99-
if (maskOIData2D.accept(i, idxNone)) {
100-
// row flagged as None:
101-
continue;
102-
}
103-
// check row flagged as Full:
104-
maskOIData2DRow = (maskOIData2D.accept(i, idxFull)) ? null : maskOIData2D;
105-
}
90+
final boolean[][] flags = oiData.getFlag();
10691

107-
final boolean[] rowFlags = (flags != null) ? flags[i] : null;
92+
IndexMask maskOIData2DRow = null;
10893

109-
// Iterate on wave channels (l):
110-
for (int l = 0; l < nWaves; l++) {
94+
// Iterate on table rows (i):
95+
for (int i = 0; i < nRows; i++) {
11196

112-
// check optional wavelength mask:
113-
if ((maskWavelength != null) && !maskWavelength.accept(l)) {
97+
// check optional data mask 1D:
98+
if ((maskOIData1D != null) && !maskOIData1D.accept(i)) {
11499
// if bit is false for this row, we hide this row
115100
continue;
116101
}
117102

118-
// check optional data mask 2D (and its Full flag):
119-
if ((maskOIData2DRow != null) && !maskOIData2DRow.accept(i, l)) {
120-
// if bit is false for this row, we hide this row
121-
continue;
103+
// check mask 2D for row None flag:
104+
if (maskOIData2D != null) {
105+
if (maskOIData2D.accept(i, idxNone)) {
106+
// row flagged as None:
107+
continue;
108+
}
109+
// check row flagged as Full:
110+
maskOIData2DRow = (maskOIData2D.accept(i, idxFull)) ? null : maskOIData2D;
122111
}
123112

124-
if ((rowFlags != null) && rowFlags[l]) {
125-
// data point is flagged so skip it:
126-
continue;
113+
final boolean[] rowFlags = (flags != null) ? flags[i] : null;
114+
115+
final double u_g = snapToGrid(ucoord[i], PREC_UV);
116+
final double v_g = snapToGrid(vcoord[i], PREC_UV);
117+
118+
// Iterate on wave channels (l):
119+
for (int l = 0; l < nWaves; l++) {
120+
121+
// check optional wavelength mask:
122+
if ((maskWavelength != null) && !maskWavelength.accept(l)) {
123+
// if bit is false for this row, we hide this row
124+
continue;
125+
}
126+
127+
// check optional data mask 2D (and its Full flag):
128+
if ((maskOIData2DRow != null) && !maskOIData2DRow.accept(i, l)) {
129+
// if bit is false for this row, we hide this row
130+
continue;
131+
}
132+
133+
if ((rowFlags != null) && rowFlags[l]) {
134+
// data point is flagged so skip it:
135+
continue;
136+
}
137+
138+
// data point is valid and not flagged:
139+
t++;
140+
final double wl_g = snapToGrid(effWaves[l], PREC_WL);
141+
142+
final double uc = u_g / wl_g; // rad-1
143+
final double vc = v_g / wl_g; // rad-1
144+
145+
// ensure (U,V) are unique (after rounding):
146+
if ((uc != 0.0) && (vc != 0.0) && uvSet.add(new UVTuple(uc, vc))) {
147+
s11 += 2.0 * (uc * uc);
148+
s22 += 2.0 * (vc * vc);
149+
s12 += 2.0 * (uc * vc);
150+
n++;
151+
}
127152
}
128-
129-
// data point is valid and not flagged:
130-
final double uc = u[i][l];
131-
final double vc = v[i][l];
132-
133-
s11 += 2.0 * (uc * uc);
134-
s22 += 2.0 * (vc * vc);
135-
s12 += 2.0 * (uc * vc);
136-
n++;
137153
}
138154
}
139155
}
156+
logger.log(Level.FINE, "t: {0}", t);
140157
logger.log(Level.FINE, "n: {0}", n);
158+
uvSet.clear();
141159

142160
if (n != 0) {
143-
// normalize by (2n + 1):
161+
// normalize by (2n + 1) to take into account symetry and (0,0) once:
144162
final double invNorm = 1.0 / (2.0 * n + 1.0);
145163

146164
s11 *= invNorm;
@@ -154,30 +172,27 @@ else if (oiData instanceof OIVis) {
154172
}
155173

156174
public static BeamInfo computeBeamInfo(final Collection<OIData> oiDatas) {
157-
int n = 0;
175+
int n = 0, t = 0;
158176
double s11 = 0.0;
159177
double s22 = 0.0;
160178
double s12 = 0.0;
161179

162-
// TODO: how to handle only OIVIS or OIT3 coords or mix (VIS/VIS2/T3) tables that may have repeated U-V coords ?
180+
final TreeSet<UVTuple> uvSet = new TreeSet<>();
181+
163182
for (OIData oiData : oiDatas) {
164-
final double[][] u; // rad-1
165-
final double[][] v; // rad-1
183+
final double[] ucoord; // m
184+
final double[] vcoord; // m
166185

167186
if (oiData instanceof OIVis2) {
168187
final OIVis2 vis2 = (OIVis2) oiData;
169-
170-
u = vis2.getSpatialUCoord();
171-
v = vis2.getSpatialVCoord();
172-
/*
173-
}
174-
else if (oiData instanceof OIVis) {
188+
ucoord = vis2.getUCoord();
189+
vcoord = vis2.getVCoord();
190+
} else if (oiData instanceof OIVis) {
175191
final OIVis vis = (OIVis) oiData;
176-
177-
u = vis.getSpatialUCoord();
178-
v = vis.getSpatialVCoord();
179-
*/
192+
ucoord = vis.getUCoord();
193+
vcoord = vis.getVCoord();
180194
} else {
195+
/** ignore OI_T3 as their UV coordinates should be redudant with OI_VIS/OI_VIS2 tables by design */
181196
continue;
182197
}
183198

@@ -186,23 +201,41 @@ else if (oiData instanceof OIVis) {
186201
final int nRows = oiData.getNbRows();
187202
final int nWaves = oiData.getNWave();
188203

189-
for (int i = 0, j; i < nRows; i++) {
190-
final boolean[] rowFlags = (flags != null) ? flags[i] : null;
204+
if (nWaves != 0) {
205+
final double[] effWaves = oiData.getOiWavelength().getEffWaveAsDouble();
206+
207+
// Iterate on table rows (i):
208+
for (int i = 0; i < nRows; i++) {
209+
final boolean[] rowFlags = (flags != null) ? flags[i] : null;
210+
211+
final double u_g = snapToGrid(ucoord[i], PREC_UV);
212+
final double v_g = snapToGrid(vcoord[i], PREC_UV);
191213

192-
for (j = 0; j < nWaves; j++) {
193-
if ((rowFlags == null) || !rowFlags[j]) {
194-
final double uc = u[i][j];
195-
final double vc = v[i][j];
214+
// Iterate on wave channels (l):
215+
for (int l = 0; l < nWaves; l++) {
196216

197-
s11 += 2.0 * (uc * uc);
198-
s22 += 2.0 * (vc * vc);
199-
s12 += 2.0 * (uc * vc);
200-
n++;
217+
if ((rowFlags == null) || !rowFlags[l]) {
218+
t++;
219+
final double wl_g = snapToGrid(effWaves[l], PREC_WL);
220+
221+
final double uc = u_g / wl_g; // rad-1
222+
final double vc = v_g / wl_g; // rad-1
223+
224+
// ensure (U,V) are unique (after rounding):
225+
if ((uc != 0.0) && (vc != 0.0) && uvSet.add(new UVTuple(uc, vc))) {
226+
s11 += 2.0 * (uc * uc);
227+
s22 += 2.0 * (vc * vc);
228+
s12 += 2.0 * (uc * vc);
229+
n++;
230+
}
231+
}
201232
}
202233
}
203234
}
204235
}
236+
logger.log(Level.FINE, "t: {0}", t);
205237
logger.log(Level.FINE, "n: {0}", n);
238+
uvSet.clear();
206239

207240
if (n != 0) {
208241
// normalize by (2n + 1):
@@ -218,6 +251,10 @@ else if (oiData instanceof OIVis) {
218251
return null;
219252
}
220253

254+
private static double snapToGrid(final double value, final double eps) {
255+
return Math.round(value / eps) * eps;
256+
}
257+
221258
private static BeamInfo computeBeam(final double[] covMatrix) {
222259

223260
if (logger.isLoggable(Level.FINE)) {
@@ -391,7 +428,10 @@ private BeamEstimator() {
391428
public static void main(String[] args) {
392429
Locale.setDefault(Locale.US);
393430
try {
394-
OIFitsFile oiFits = OIFitsLoader.loadOIFits("/home/bourgesl/Documents/vlti2023/data/vlti/jmmctools/ImageReconstruction/data/RCar/R_CAR_all.fits");
431+
final OIFitsFile oiFits = OIFitsLoader.loadOIFits(
432+
// "/home/bourgesl/Documents/vlti2023/data/vlti/jmmctools/ImageReconstruction/data/RCar/R_CAR_all.fits"
433+
"/home/bourgesl/ASPRO2/oifits/Aspro2_HIP1234_VLTI_MATISSE_LM_2.86542-4.18239-62ch_D0-G2-J3-K0_2023-10-29.fits"
434+
);
395435

396436
final BeamInfo beamInfo = computeBeamInfo(oiFits.getOiDataList());
397437

@@ -456,4 +496,55 @@ public void paint(Graphics g) {
456496
frame.setVisible(true);
457497
}
458498

499+
final static class UVTuple implements Comparable<UVTuple> {
500+
501+
final double u;
502+
final double v;
503+
504+
UVTuple(double u, double v) {
505+
this.u = u;
506+
this.v = v;
507+
}
508+
509+
@Override
510+
public int hashCode() {
511+
int hash = 7;
512+
hash = 11 * hash + (int) (Double.doubleToLongBits(this.u) ^ (Double.doubleToLongBits(this.u) >>> 32));
513+
hash = 11 * hash + (int) (Double.doubleToLongBits(this.v) ^ (Double.doubleToLongBits(this.v) >>> 32));
514+
return hash;
515+
}
516+
517+
@Override
518+
public boolean equals(Object obj) {
519+
if (this == obj) {
520+
return true;
521+
}
522+
if (obj == null) {
523+
return false;
524+
}
525+
if (getClass() != obj.getClass()) {
526+
return false;
527+
}
528+
final UVTuple other = (UVTuple) obj;
529+
if (Double.doubleToLongBits(this.u) != Double.doubleToLongBits(other.u)) {
530+
return false;
531+
}
532+
return Double.doubleToLongBits(this.v) == Double.doubleToLongBits(other.v);
533+
}
534+
535+
@Override
536+
public int compareTo(final UVTuple o) {
537+
int res = Double.compare(u, o.u);
538+
if (res == 0) {
539+
res = Double.compare(v, o.v);
540+
}
541+
return res;
542+
}
543+
544+
@Override
545+
public String toString() {
546+
return "(" + "u=" + u + ", v=" + v + ')';
547+
}
548+
549+
}
459550
}

0 commit comments

Comments
 (0)