Source
309
309
Bool MultiTermMatrixCleaner::getinvhessian(Matrix<Double> & invhessian)
310
310
{
311
311
invhessian.resize((invMatA_p[0]).shape());
312
312
invhessian = (invMatA_p[0]);
313
313
return true;
314
314
}
315
315
316
316
317
317
/* Do the deconvolution */
318
318
Int MultiTermMatrixCleaner::mtclean(Int maxniter, Float stopfraction, Float inputgain, Float userthreshold)
319
-
{
319
+
{
320
320
LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
321
321
if(adbg)os << "SOLVER for Multi-Frequency Synthesis deconvolution" << LogIO::POST;
322
322
323
323
/* Initialize some variables */
324
324
325
325
maxniter_p = maxniter;
326
326
stopfraction_p= stopfraction;
327
327
inputgain_p=inputgain;
328
328
userthreshold_p=userthreshold;
329
329
357
357
/* Check if peak RHS is already less than the threshold */
358
358
/* This step computes the flux limit also - when called here... */
359
359
convergedflag = checkConvergence(criterion,fluxlimit,loopgain);
360
360
if(convergedflag==1) return 0;
361
361
362
362
/* Compute the flux limits that determine the depth of the minor cycles. */
363
363
// computeFluxLimit(fluxlimit,thresh);
364
364
365
365
/********************** START MINOR CYCLE ITERATIONS ***********************/
366
366
//os << "Doing deconvolution iterations..." << LogIO::POST;
367
+
367
368
for(itercount_p=0;itercount_p<maxniter_p;itercount_p++)
368
369
{
369
370
globalmaxval_p=-1e+10;
370
371
371
372
// if itercount_p==0, set the blc/trc to the full image size.
372
373
// For later iterations, it will get updated according to globalmaxpos_p and psfsupport_p
373
374
// when buildImagePatches() is called.
374
375
if(itercount_p==0)
375
376
{
376
377
blc_p = IPosition(2,0,0);
392
393
// Find max across all pixels.
393
394
chooseComponent(ntaylor, scale,criterion,blc,trc);
394
395
}
395
396
}//end pragma omp
396
397
397
398
/* Find the best component over all scales */
398
399
for(Int scale=0;scale<nscales_p;scale++)
399
400
{
400
401
if((maxScaleVal_p[scale]*scaleBias_p[scale]) > globalmaxval_p)
401
402
{
402
-
globalmaxval_p = maxScaleVal_p[scale];
403
+
globalmaxval_p = maxScaleVal_p[scale]*scaleBias_p[scale];
403
404
globalmaxpos_p = maxScalePos_p[scale];
404
-
maxscaleindex_p = scale;
405
+
maxscaleindex_p = scale;
405
406
}
406
407
}
407
408
408
409
/* Update the image and psf patch sizes according to the
409
410
current globalmaxval and psfsupport.
410
411
This patch is over which the next-iteration's solveMatrixEqn is computed */
411
412
buildImagePatches();
412
413
413
414
414
415
/* Update the current solution by this chosen step */