Source
xxxxxxxxxx
2653
2653
2654
2654
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2655
2655
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2656
2656
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2657
2657
/////////////////// Utility Functions to gather statistics on the imagestore.
2658
2658
2659
2659
Float SIImageStore::getPeakResidual()
2660
2660
{
2661
2661
LogIO os( LogOrigin("SIImageStore","getPeakResidual",WHERE) );
2662
2662
2663
-
//Lattice<Bool> pixelmask = *residual()->pixelMask();
2664
-
//LatticeExprNode pres( max(abs( *residual() ) ));
2665
-
//Lattice<Bool> pixelmask = *residual()->pixelMask();
2666
-
os<<"hasPixelMask="<<residual()->hasPixelMask()<<LogIO::POST;
2667
-
os<<"ntrue="<< ntrue(residual()->getMask()) << LogIO::POST;
2668
2663
ArrayLattice<Bool> pixelmask(residual()->getMask());
2669
2664
LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
2670
2665
//LatticeExprNode pres( max(abs( *residual() ) ));
2671
2666
LatticeExprNode pres( max(abs(resd) ));
2672
2667
Float maxresidual = pres.getFloat();
2673
2668
2674
-
// Float maxresidual = max( residual()->get() );
2675
-
//Record* regionPtr=0;
2676
-
//String LELmask("");
2677
-
//Record thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
2678
-
//Array<Double> maxs, mins;
2679
-
//Double maxofabsmaxs = max(abs(maxs));
2680
-
//os<<"mxofabsmaxs ===="<<maxofabsmaxs<<LogIO::POST;
2681
-
//Double maxofabsmins = max(abs(mins));
2682
-
//Float maxresidual = (Float) max(maxofabsmaxs,maxofabsmins);
2683
2669
2684
2670
return maxresidual;
2685
2671
}
2686
2672
2687
2673
Float SIImageStore::getPeakResidualWithinMask()
2688
2674
{
2689
2675
LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
2690
2676
Float minresmask, maxresmask, minres, maxres;
2691
2677
//findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2692
2678
2693
-
os<<"getPeakResWithinMask hasPixelMask="<<residual()->hasPixelMask()<<LogIO::POST;
2694
-
os<<"ntrue="<< ntrue(residual()->getMask()) <<LogIO::POST;
2695
2679
ArrayLattice<Bool> pixelmask(residual()->getMask());
2696
2680
findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2697
2681
2698
2682
//return maxresmask;
2699
2683
return max( abs(maxresmask), abs(minresmask) );
2700
2684
}
2701
2685
2702
2686
// Calculate the total model flux
2703
2687
Float SIImageStore::getModelFlux(uInt term)
2704
2688
{
2869
2853
Float minresmask=0, maxresmask=0, minres=0, maxres=0;
2870
2854
// findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2871
2855
ArrayLattice<Bool> pixelmask(residual()->getMask());
2872
2856
if(hasMask())
2873
2857
{
2874
2858
findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2875
2859
}
2876
2860
else
2877
2861
{
2878
2862
//LatticeExprNode pres( max( *residual() ) );
2879
-
LatticeExprNode pres( max( iif(!pixelmask,*residual(),0) ) );
2863
+
LatticeExprNode pres( max( iif(pixelmask,*residual(),0) ) );
2880
2864
maxres = pres.getFloat();
2881
2865
//LatticeExprNode pres2( min( *residual() ) );
2882
-
LatticeExprNode pres2( min( iif(!pixelmask,*residual(),0) ) );
2866
+
LatticeExprNode pres2( min( iif(pixelmask,*residual(),0) ) );
2883
2867
minres = pres2.getFloat();
2884
2868
}
2885
2869
2886
2870
os << "[" << itsImageName << "]" ;
2887
2871
os << " Peak residual (max,min) " ;
2888
2872
if( minresmask!=0.0 || maxresmask!=0.0 )
2889
2873
{ os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
2890
2874
os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
2891
2875
2892
2876
os << "[" << itsImageName << "] Total Model Flux : " << getModelFlux() << LogIO::POST;
2893
2877
2894
2878
2895
2879
Record* regionPtr=0;
2896
2880
String LELmask("");
2897
2881
Record thestats = SDMaskHandler::calcImageStatistics(*residual(), LELmask, regionPtr, True);
2898
2882
Array<Double> maxs, mins;
2899
2883
thestats.get(RecordFieldId("max"), maxs);
2900
2884
thestats.get(RecordFieldId("min"), mins);
2901
-
os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
2902
-
os << LogIO::DEBUG1 << "Min : " << mins << LogIO::POST;
2885
+
//os << LogIO::DEBUG1 << "Max : " << maxs << LogIO::POST;
2886
+
//os << LogIO::DEBUG1 << "Min : " << mins << LogIO::POST;
2903
2887
2904
2888
}
2905
2889
2906
2890
// Calculate the total model flux
2907
2891
Float SIImageStore::getMaskSum()
2908
2892
{
2909
2893
LogIO os( LogOrigin("SIImageStore","getMaskSum",WHERE) );
2910
2894
2911
2895
LatticeExprNode msum( sum( *mask() ) );
2912
2896
Float masksum = msum.getFloat();
2917
2901
}
2918
2902
2919
2903
Bool SIImageStore::findMinMaxLattice(const Lattice<Float>& lattice,
2920
2904
const Lattice<Float>& mask,
2921
2905
const Lattice<Bool>& pixelmask,
2922
2906
Float& maxAbs, Float& maxAbsMask,
2923
2907
Float& minAbs, Float& minAbsMask )
2924
2908
{
2925
2909
2926
2910
//FOR DEGUG
2927
-
LogIO os( LogOrigin("SIImageStore","findMinMaxLattice",WHERE) );
2911
+
//LogIO os( LogOrigin("SIImageStore","findMinMaxLattice",WHERE) );
2928
2912
2929
2913
maxAbs=0.0;maxAbsMask=0.0;
2930
2914
minAbs=1e+10;minAbsMask=1e+10;
2931
2915
2932
2916
const IPosition tileShape = lattice.niceCursorShape();
2933
2917
TiledLineStepper ls(lattice.shape(), tileShape, 0);
2934
2918
{
2935
2919
RO_LatticeIterator<Float> li(lattice, ls);
2936
2920
RO_LatticeIterator<Float> mi(mask, ls);
2937
2921
RO_LatticeIterator<Bool> pmi(pixelmask, ls);
2938
2922
for(li.reset(),mi.reset(),pmi.reset();!li.atEnd();li++, mi++, pmi++) {
2939
2923
IPosition posMax=li.position();
2940
2924
IPosition posMin=li.position();
2941
2925
IPosition posMaxMask=li.position();
2942
2926
IPosition posMinMask=li.position();
2943
2927
Float maxVal=0.0;
2944
2928
Float minVal=0.0;
2945
2929
Float maxValMask=0.0;
2946
2930
Float minValMask=0.0;
2947
2931
2948
-
os<<"nelements for lattice="<<(li.cursor()).nelements()<<LogIO::POST;
2949
-
os<<"nelements for pmi="<<(pmi.cursor()).nelements()<<LogIO::POST;
2950
-
os<<"ntrue for pmi="<<ntrue(pmi.cursor())<<LogIO::POST;
2951
2932
2952
2933
// skip if lattice chunk is masked entirely.
2953
2934
if(ntrue(pmi.cursor()) > 0 ) {
2954
2935
MaskedArray<Float> marr(li.cursor(), pmi.cursor());
2955
2936
MaskedArray<Float> marrinmask(li.cursor() * mi.cursor(), pmi.cursor());
2956
2937
//minMax( minVal, maxVal, posMin, posMax, li.cursor() );
2957
-
//minMax( minVal, maxVal, posMin, posMax, marr );
2958
-
os<<"do minMax with marrinmask.nelements()="<<marrinmask.nelements()<<LogIO::POST;
2938
+
minMax( minVal, maxVal, posMin, posMax, marr );
2959
2939
//minMaxMasked(minValMask, maxValMask, posMin, posMax, li.cursor(), mi.cursor());
2960
2940
minMax(minValMask, maxValMask, posMin, posMax, marrinmask);
2961
2941
2962
-
os<<"DONE minMax"<<LogIO::POST;
2942
+
//os<<"DONE minMax"<<LogIO::POST;
2963
2943
if( (maxVal) > (maxAbs) ) maxAbs = maxVal;
2964
2944
if( (maxValMask) > (maxAbsMask) ) maxAbsMask = maxValMask;
2965
2945
2966
2946
if( (minVal) < (minAbs) ) minAbs = minVal;
2967
2947
if( (minValMask) < (minAbsMask) ) minAbsMask = minValMask;
2968
2948
}
2969
2949
2970
2950
}
2971
2951
}
2972
2952