Commits

Renaud Miel authored e2744a57573
CAS-10682 Removed slicing loop in SDGrid::makeImage

Restuctured and cleaned up code
No tags

casa5/code/synthesis/MeasurementComponents/SDGrid.cc

Modified
1240 1240 &convSampling,
1241 1241 convFunc.getStorage(del),
1242 1242 chanMap.getStorage(del),
1243 1243 polMap.getStorage(del));
1244 1244
1245 1245 data.putStorage(datStorage, isCopy);
1246 1246 }
1247 1247 interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1248 1248 }
1249 1249
1250 -
1251 -
1252 - // Make a plain straightforward honest-to-FSM image. This returns
1253 - // a complex image, without conversion to Stokes. The representation
1254 - // is that required for the visibilities.
1255 - //----------------------------------------------------------------------
1256 - void SDGrid::makeImage(FTMachine::Type type,
1257 - ROVisibilityIterator& vi,
1258 - ImageInterface<Complex>& theImage,
1259 - Matrix<Float>& weight) {
1260 -
1250 +// Make a plain straightforward honest-to-FSM image. This returns
1251 +// a complex image, without conversion to Stokes. The representation
1252 +// is that required for the visibilities.
1253 +//----------------------------------------------------------------------
1254 +void SDGrid::makeImage(FTMachine::Type inType,
1255 + ROVisibilityIterator& vi,
1256 + ImageInterface<Complex>& theImage,
1257 + Matrix<Float>& weight) {
1261 1258
1262 1259 logIO() << LogOrigin("FTMachine", "makeImage0") << LogIO::NORMAL;
1263 1260
1264 - // Loop over all visibilities and pixels
1261 + // Attach buffer to iterator
1265 1262 VisBuffer vb(vi);
1266 1263
1267 - // Initialize put (i.e. transform to Sky) for this model
1264 + // Set CStokesRep
1268 1265 vi.origin();
1269 -
1270 - if(vb.polFrame()==MSIter::Linear) {
1271 - StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::LINEAR);
1272 - }
1273 - else {
1274 - StokesImageUtil::changeCStokesRep(theImage, StokesImageUtil::CIRCULAR);
1275 - }
1276 - Bool useCorrected= !(vi.msColumns().correctedData().isNull());
1277 - if((type==FTMachine::CORRECTED) && (!useCorrected))
1278 - type=FTMachine::OBSERVED;
1279 - Bool normalize=true;
1280 - if(type==FTMachine::COVERAGE)
1281 - normalize=false;
1282 -
1283 - Int Nx=theImage.shape()(0);
1284 - Int Ny=theImage.shape()(1);
1285 - Int Npol=theImage.shape()(2);
1286 - Int Nchan=theImage.shape()(3);
1287 - Double memtot=Double(HostInfo::memoryTotal(true))*1024.0; // return in kB
1288 - Int nchanInMem=Int(memtot/2.0/8.0/Double(Nx*Ny*Npol));
1289 - Int nloop=nchanInMem >= Nchan ? 1 : Nchan/nchanInMem+1;
1290 - ImageInterface<Complex> *imCopy=NULL;
1291 - IPosition blc(theImage.shape());
1292 - IPosition trc(theImage.shape());
1293 - blc-=blc; //set all values to 0
1294 - trc=theImage.shape();
1295 - trc-=1; // set trc to image size -1
1296 - if(nloop==1) {
1297 - imCopy=& theImage;
1298 - nchanInMem=Nchan;
1299 - }
1300 - else
1301 - logIO() << "Not enough memory to image in one go \n will process the image in "
1302 - << nloop
1303 - << " sections "
1304 - << LogIO::POST;
1305 -
1266 + auto cStokesRep = (vb.polFrame() == MSIter::Linear) ?
1267 + StokesImageUtil::LINEAR : StokesImageUtil::CIRCULAR;
1268 + StokesImageUtil::changeCStokesRep(theImage, cStokesRep);
1269 +
1270 + // Weights
1271 + auto imageShape = theImage.shape();
1272 + auto Npol = imageShape(2);
1273 + auto Nchan = imageShape(3);
1306 1274 weight.resize(Npol, Nchan);
1307 - Matrix<Float> wgtcopy(Npol, Nchan);
1308 -
1309 - Bool isWgtZero=true;
1310 - for (Int k=0; k < nloop; ++k){
1311 - Int bchan=k*nchanInMem;
1312 - Int echan=(k+1)*nchanInMem < Nchan ? (k+1)*nchanInMem-1 : Nchan-1;
1313 -
1314 - if(nloop > 1) {
1315 - blc[3]=bchan;
1316 - trc[3]=echan;
1317 - Slicer sl(blc, trc, Slicer::endIsLast);
1318 - imCopy=new SubImage<Complex>(theImage, sl, true);
1319 - wgtcopy.resize(npol, echan-bchan+1);
1320 - }
1321 - vi.originChunks();
1322 - vi.origin();
1323 - initializeToSky(*imCopy,wgtcopy,vb);
1324 1275
1276 + initializeToSky(theImage,weight,vb);
1277 + // Loop over the visibilities, putting VisBuffers
1278 + for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1279 + for (vi.origin(); vi.more(); vi++) {
1280 + FTMachine::Type actualType;
1281 + Bool doPSF;
1282 + getParamsForFTMachineType(vi, inType, doPSF, actualType);
1283 + setupVisBufferForFTMachineType(actualType,vb);
1284 + constexpr Int allVbRows = -1;
1285 + put(vb, allVbRows, doPSF, actualType);
1286 + }
1287 + }
1288 + finalizeToSky();
1325 1289
1326 - // for minmax clipping
1327 - logIO() << LogOrigin("SDGrid", "makeImage", WHERE) << LogIO::DEBUGGING
1328 - << "doclip_ = " << (clipminmax_ ? "TRUE" : "FALSE") << " (" << clipminmax_ << ")" << LogIO::POST;
1329 - if (clipminmax_) {
1330 - logIO() << LogOrigin("SDGRID", "makeImage", WHERE)
1331 - << LogIO::DEBUGGING << "use ggridsd2 for imaging" << LogIO::POST;
1332 - }
1290 + // Normalize by dividing out weights, etc.
1291 + auto doNormalize = (inType == FTMachine::COVERAGE) ? false : true;
1292 + getImage(weight, doNormalize);
1333 1293
1334 - // Loop over the visibilities, putting VisBuffers
1335 - for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
1336 - for (vi.origin(); vi.more(); vi++) {
1337 -
1338 - switch(type) {
1339 - case FTMachine::RESIDUAL:
1340 - vb.visCube()=vb.correctedVisCube();
1341 - vb.visCube()-=vb.modelVisCube();
1342 - put(vb, -1, false);
1343 - break;
1344 - case FTMachine::MODEL:
1345 - put(vb, -1, false, FTMachine::MODEL);
1346 - break;
1347 - case FTMachine::CORRECTED:
1348 - put(vb, -1, false, FTMachine::CORRECTED);
1349 - break;
1350 - case FTMachine::PSF:
1351 - vb.visCube()=Complex(1.0,0.0);
1352 - put(vb, -1, true, FTMachine::PSF);
1353 - break;
1354 - case FTMachine::COVERAGE:
1355 - vb.visCube()=Complex(1.0);
1356 - put(vb, -1, true, FTMachine::COVERAGE);
1357 - break;
1358 - case FTMachine::OBSERVED:
1359 - default:
1360 - put(vb, -1, false, FTMachine::OBSERVED);
1361 - break;
1362 - }
1363 - }
1364 - }
1365 - finalizeToSky();
1366 - // Normalize by dividing out weights, etc.
1367 - getImage(wgtcopy, normalize);
1368 - if(max(wgtcopy)==0.0){
1369 - if(nloop > 1)
1370 - logIO() << LogIO::WARN
1371 - << "No useful data in SDGrid: weights all zero for image slice " << k
1372 - << LogIO::POST;
1373 - }
1374 - else
1375 - isWgtZero=false;
1376 -
1377 - weight(Slice(0, Npol), Slice(bchan, echan-bchan+1))=wgtcopy;
1378 - if(nloop >1) delete imCopy;
1379 - }//loop k
1380 - if(isWgtZero)
1381 - logIO() << LogIO::SEVERE
1382 - << "No useful data in SDGrid: weights all zero"
1383 - << LogIO::POST;
1384 - }
1294 + // Warning message
1295 + if (allEQ(weight,0.0f)) {
1296 + logIO() << LogIO::SEVERE
1297 + << "No useful data in SDGrid: all weights are zero"
1298 + << LogIO::POST;
1299 + }
1300 +}
1385 1301
1302 +void SDGrid::getParamsForFTMachineType(const ROVisibilityIterator& vi, FTMachine::Type in_type,
1303 + casacore::Bool& out_dopsf, FTMachine::Type& out_type) const {
1304 +
1305 + // Tune input type of FTMachine
1306 + auto haveCorrectedData = not (vi.msColumns().correctedData().isNull());
1307 + auto tunedType =
1308 + ((in_type == FTMachine::CORRECTED) && (not haveCorrectedData)) ?
1309 + FTMachine::OBSERVED : in_type;
1310 +
1311 + // Compute output parameters
1312 + switch(tunedType) {
1313 + case FTMachine::RESIDUAL:
1314 + case FTMachine::MODEL:
1315 + case FTMachine::CORRECTED:
1316 + out_dopsf = false;
1317 + out_type = tunedType;
1318 + break;
1319 + case FTMachine::PSF:
1320 + case FTMachine::COVERAGE:
1321 + out_dopsf = true;
1322 + out_type = tunedType;
1323 + break;
1324 + default:
1325 + out_dopsf = false;
1326 + out_type = FTMachine::OBSERVED;
1327 + break;
1328 + }
1329 +}
1386 1330
1331 +void SDGrid::setupVisBufferForFTMachineType(FTMachine::Type type, VisBuffer& vb) const {
1332 + switch(type) {
1333 + case FTMachine::RESIDUAL:
1334 + vb.visCube() = vb.correctedVisCube();
1335 + vb.visCube() -= vb.modelVisCube();
1336 + break;
1337 + case FTMachine::PSF:
1338 + vb.visCube()=Complex(1.0,0.0);
1339 + break;
1340 + case FTMachine::COVERAGE:
1341 + vb.visCube()=Complex(1.0);
1342 + break;
1343 + default:
1344 + break;
1345 + }
1346 +}
1387 1347
1388 1348 // Finalize : optionally normalize by weight image
1389 1349 ImageInterface<Complex>& SDGrid::getImage(Matrix<Float>& weights,
1390 1350 Bool normalize)
1391 1351 {
1392 1352 AlwaysAssert(lattice, AipsError);
1393 1353 AlwaysAssert(image, AipsError);
1394 1354
1395 1355 logIO() << LogOrigin("SDGrid", "getImage") << LogIO::NORMAL;
1396 1356
1552 1512 }
1553 1513 }
1554 1514 // No match!
1555 1515 return -1;
1556 1516 }
1557 1517
1558 1518 Bool SDGrid::getXYPos(const VisBuffer& vb, Int row) {
1559 1519
1560 1520 Bool dointerp;
1561 1521 Bool nullPointingTable = false;
1522 + // const MSPointingColumns& act_mspc =
1523 + // interpolationRefFrame == InterpolationRefFrame::POINTING ?
1524 + // vb.msColumns().pointing() : convertedPointing_p ;
1562 1525 const MSPointingColumns& act_mspc = vb.msColumns().pointing();
1563 1526 nullPointingTable = (act_mspc.nrow() < 1);
1564 1527 Int pointIndex = -1;
1565 1528 if (!nullPointingTable) {
1566 1529 ///if(vb.newMS()) vb.newMS does not work well using msid
1567 1530 if (vb.msId() != msId_p) {
1568 1531 lastIndex_p = 0;
1569 1532 if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1570 1533 lastIndexPerAnt_p.resize(vb.numberAnt());
1571 1534 }

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut