Commits
Takeshi Nakazato authored and Ville Suoranta committed 73c121f99c0 Merge
2153 2153 | void SDGrid::pickWeights(const VisBuffer& vb, Matrix<Float>& weight){ |
2154 2154 | |
2155 2155 | StartStop trigger(cPickWeights); |
2156 2156 | |
2157 2157 | //break reference |
2158 2158 | weight.resize(); |
2159 2159 | |
2160 2160 | if (useImagingWeight_p) { |
2161 2161 | weight.reference(vb.imagingWeight()); |
2162 2162 | } else { |
2163 - | const Cube<Float> weightspec(vb.weightSpectrum()); |
2163 + | const Cube<Float> weightSpec(vb.weightSpectrum()); |
2164 2164 | weight.resize(vb.nChannel(), vb.nRow()); |
2165 2165 | |
2166 - | if (weightspec.nelements() == 0) { |
2167 - | for (Int k = 0; k < vb.nRow(); ++k) { |
2168 - | //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl; |
2169 - | weight.column(k).set(vb.weight()(k)); |
2166 + | // CAS-9957 correct weight propagation from linear/circular correlations to Stokes I |
2167 + | const auto toStokesWeight = [](float weight0, float weight1) { |
2168 + | const auto denominator = weight0 + weight1; |
2169 + | const auto numerator = weight0 * weight1; |
2170 + | constexpr float fmin = std::numeric_limits<float>::min(); |
2171 + | return abs(denominator) < fmin ? 0.0f : 4.0f * numerator / denominator; |
2172 + | }; |
2173 + | |
2174 + | if (weightSpec.nelements() == 0) { |
2175 + | const auto &weightMat = vb.weightMat(); |
2176 + | const ssize_t npol = weightMat.shape()(0); |
2177 + | if (npol == 1) { |
2178 + | const auto weight0 = weightMat.row(0); |
2179 + | for (int k = 0; k < vb.nRow(); ++k) { |
2180 + | weight.column(k).set(weight0(k)); |
2181 + | } |
2182 + | } else if (npol == 2) { |
2183 + | const auto weight0 = weightMat.row(0); |
2184 + | const auto weight1 = weightMat.row(1); |
2185 + | for (int k = 0; k < vb.nRow(); ++k) { |
2186 + | //cerr << "nrow " << vb.nRow() << " " << weight.shape() << " " << weight.column(k).shape() << endl; |
2187 + | weight.column(k).set(toStokesWeight(weight0(k), weight1(k))); |
2188 + | } |
2189 + | } else { |
2190 + | // It seems current code doesn't support 4 pol case. So, give up |
2191 + | // processing such data to avoid producing unintended result |
2192 + | throw AipsError("Imaging full-Stokes data (npol=4) is not supported."); |
2170 2193 | } |
2171 2194 | } else { |
2172 - | Int npol = weightspec.shape()(0); |
2195 + | const ssize_t npol = weightSpec.shape()(0); |
2173 2196 | if (npol == 1) { |
2174 - | for (Int k = 0; k < vb.nRow(); ++k) { |
2197 + | weight = weightSpec.yzPlane(0); |
2198 + | } else if (npol == 2) { |
2199 + | const auto weight0 = weightSpec.yzPlane(0); |
2200 + | const auto weight1 = weightSpec.yzPlane(1); |
2201 + | for (int k = 0; k < vb.nRow(); ++k) { |
2175 2202 | for (int chan = 0; chan < vb.nChannel(); ++chan) { |
2176 - | weight(chan, k)=weightspec(0, chan, k); |
2203 + | weight(chan, k) = toStokesWeight(weight0(chan, k), weight1(chan, k)); |
2177 2204 | } |
2178 2205 | } |
2179 2206 | } else { |
2180 - | for (Int k = 0; k < vb.nRow(); ++k) { |
2181 - | for (int chan = 0; chan < vb.nChannel(); ++chan) { |
2182 - | weight(chan, k) = (weightspec(0, chan, k) + weightspec((npol-1), chan, k))/2.0f; |
2183 - | } |
2184 - | } |
2207 + | // It seems current code doesn't support 4 pol case. So, give up |
2208 + | // processing such data to avoid producing unintended result |
2209 + | throw AipsError("Imaging full-Stokes data (npol=4) is not supported."); |
2185 2210 | } |
2186 2211 | } |
2187 2212 | } |
2188 2213 | } |
2189 2214 | |
2190 2215 | void SDGrid::clipMinMax() { |
2191 2216 | if (clipminmax_) { |
2192 2217 | Bool gmin_b, gmax_b, wmin_b, wmax_b, np_b; |
2193 2218 | const auto *gmin_p = gmin_.getStorage(gmin_b); |
2194 2219 | const auto *gmax_p = gmax_.getStorage(gmax_b); |