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
-
LatticeExprNode pres( max(abs( *residual() ) ));
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
+
ArrayLattice<Bool> pixelmask(residual()->getMask());
2669
+
LatticeExpr<Float> resd(iif(pixelmask,abs(*residual()),0));
2670
+
//LatticeExprNode pres( max(abs( *residual() ) ));
2671
+
LatticeExprNode pres( max(abs(resd) ));
2664
2672
Float maxresidual = pres.getFloat();
2665
2673
2666
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);
2667
2683
2668
2684
return maxresidual;
2669
2685
}
2670
2686
2671
2687
Float SIImageStore::getPeakResidualWithinMask()
2672
2688
{
2673
2689
LogIO os( LogOrigin("SIImageStore","getPeakResidualWithinMask",WHERE) );
2674
2690
Float minresmask, maxresmask, minres, maxres;
2675
2691
//findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2676
2692
2677
-
findMinMaxLattice(*residual(), *mask() , maxres,maxresmask, minres, minresmask);
2693
+
os<<"getPeakResWithinMask hasPixelMask="<<residual()->hasPixelMask()<<LogIO::POST;
2694
+
os<<"ntrue="<< ntrue(residual()->getMask()) <<LogIO::POST;
2695
+
ArrayLattice<Bool> pixelmask(residual()->getMask());
2696
+
findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2678
2697
2679
2698
//return maxresmask;
2680
2699
return max( abs(maxresmask), abs(minresmask) );
2681
2700
}
2682
2701
2683
2702
// Calculate the total model flux
2684
2703
Float SIImageStore::getModelFlux(uInt term)
2685
2704
{
2686
2705
// LogIO os( LogOrigin("SIImageStore","getModelFlux",WHERE) );
2687
2706
2842
2861
//return mdns+mads*1.4826;
2843
2862
// this is the old noise calc
2844
2863
return mads*1.4826;
2845
2864
}
2846
2865
2847
2866
void SIImageStore::printImageStats()
2848
2867
{
2849
2868
LogIO os( LogOrigin("SIImageStore","printImageStats",WHERE) );
2850
2869
Float minresmask=0, maxresmask=0, minres=0, maxres=0;
2851
2870
// findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
2871
+
ArrayLattice<Bool> pixelmask(residual()->getMask());
2852
2872
if(hasMask())
2853
2873
{
2854
-
findMinMaxLattice(*residual(), *mask() , maxres,maxresmask, minres, minresmask);
2874
+
findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
2855
2875
}
2856
2876
else
2857
2877
{
2858
-
LatticeExprNode pres( max( *residual() ) );
2878
+
//LatticeExprNode pres( max( *residual() ) );
2879
+
LatticeExprNode pres( max( iif(!pixelmask,*residual(),0) ) );
2859
2880
maxres = pres.getFloat();
2860
-
LatticeExprNode pres2( min( *residual() ) );
2881
+
//LatticeExprNode pres2( min( *residual() ) );
2882
+
LatticeExprNode pres2( min( iif(!pixelmask,*residual(),0) ) );
2861
2883
minres = pres2.getFloat();
2862
2884
}
2863
2885
2864
2886
os << "[" << itsImageName << "]" ;
2865
2887
os << " Peak residual (max,min) " ;
2866
2888
if( minresmask!=0.0 || maxresmask!=0.0 )
2867
2889
{ os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
2868
2890
os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
2869
2891
2870
2892
os << "[" << itsImageName << "] Total Model Flux : " << getModelFlux() << LogIO::POST;
2889
2911
LatticeExprNode msum( sum( *mask() ) );
2890
2912
Float masksum = msum.getFloat();
2891
2913
2892
2914
// Float masksum = sum( mask()->get() );
2893
2915
2894
2916
return masksum;
2895
2917
}
2896
2918
2897
2919
Bool SIImageStore::findMinMaxLattice(const Lattice<Float>& lattice,
2898
2920
const Lattice<Float>& mask,
2921
+
const Lattice<Bool>& pixelmask,
2899
2922
Float& maxAbs, Float& maxAbsMask,
2900
2923
Float& minAbs, Float& minAbsMask )
2901
2924
{
2902
2925
2926
+
//FOR DEGUG
2927
+
LogIO os( LogOrigin("SIImageStore","findMinMaxLattice",WHERE) );
2928
+
2903
2929
maxAbs=0.0;maxAbsMask=0.0;
2904
2930
minAbs=1e+10;minAbsMask=1e+10;
2905
2931
2906
2932
const IPosition tileShape = lattice.niceCursorShape();
2907
2933
TiledLineStepper ls(lattice.shape(), tileShape, 0);
2908
2934
{
2909
2935
RO_LatticeIterator<Float> li(lattice, ls);
2910
2936
RO_LatticeIterator<Float> mi(mask, ls);
2911
-
for(li.reset(),mi.reset();!li.atEnd();li++, mi++) {
2937
+
RO_LatticeIterator<Bool> pmi(pixelmask, ls);
2938
+
for(li.reset(),mi.reset(),pmi.reset();!li.atEnd();li++, mi++, pmi++) {
2912
2939
IPosition posMax=li.position();
2913
2940
IPosition posMin=li.position();
2914
2941
IPosition posMaxMask=li.position();
2915
2942
IPosition posMinMask=li.position();
2916
2943
Float maxVal=0.0;
2917
2944
Float minVal=0.0;
2918
2945
Float maxValMask=0.0;
2919
2946
Float minValMask=0.0;
2920
-
2921
-
minMaxMasked(minValMask, maxValMask, posMin, posMax, li.cursor(), mi.cursor());
2922
2947
2923
-
minMax( minVal, maxVal, posMin, posMax, li.cursor() );
2924
-
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
+
2952
+
// skip if lattice chunk is masked entirely.
2953
+
if(ntrue(pmi.cursor()) > 0 ) {
2954
+
MaskedArray<Float> marr(li.cursor(), pmi.cursor());
2955
+
MaskedArray<Float> marrinmask(li.cursor() * mi.cursor(), pmi.cursor());
2956
+
//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;
2959
+
//minMaxMasked(minValMask, maxValMask, posMin, posMax, li.cursor(), mi.cursor());
2960
+
minMax(minValMask, maxValMask, posMin, posMax, marrinmask);
2961
+
2962
+
os<<"DONE minMax"<<LogIO::POST;
2925
2963
if( (maxVal) > (maxAbs) ) maxAbs = maxVal;
2926
2964
if( (maxValMask) > (maxAbsMask) ) maxAbsMask = maxValMask;
2927
2965
2928
2966
if( (minVal) < (minAbs) ) minAbs = minVal;
2929
2967
if( (minValMask) < (minAbsMask) ) minAbsMask = minValMask;
2968
+
}
2930
2969
2931
2970
}
2932
2971
}
2933
2972
2934
2973
return True;
2935
2974
2936
2975
2937
2976
}
2938
2977
2939
2978
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////