Initial work for Image Guiding algorithm by Bob Majewski. All guide images...

Initial work for Image Guiding algorithm by Bob Majewski. All guide images must be explicitly converted to floats and then partitioned into NxN regions before sending them to the algorithm

kstars/CMakeLists.txt

@@ -202,6 +202,7 @@ if (INDI_FOUND)
                        ekos/guide/internalguide/matr.cpp
                        #ekos/guide/internalguide/rcalibration.cpp
                        ekos/guide/internalguide/vect.cpp
+                       ekos/guide/internalguide/imageautoguiding.cpp
                        # External Guide
                        ekos/guide/externalguide/phd2.cpp
                        ekos/guide/externalguide/linguider.cpp

kstars/ekos/guide/guide.cpp

@@ -1441,6 +1441,8 @@ bool Guide::setGuiderType(int type)
 
             guider = internalGuider;
 
+            internalGuider->setRegionAxis(opsGuide->kcfg_GuideRegionAxis->currentText().toInt());
+
             calibrateB->setEnabled(true);
             guideB->setEnabled(false);
             captureB->setEnabled(true);
@@ -1985,8 +1987,13 @@ void Guide::buildOperationStack(GuideState operation)
                 if (subFramed == false && Options::guideSubframeEnabled() && Options::guideAutoStarEnabled())
                     operationStack.push(GUIDE_CAPTURE);
 
-                operationStack.push(GUIDE_SUBFRAME);
-                operationStack.push(GUIDE_STAR_SELECT);
+                // Do not subframe and auto-select star on Image Guiding mode
+                if (Options::imageGuidingEnabled() == false)
+                {
+                    operationStack.push(GUIDE_SUBFRAME);
+                    operationStack.push(GUIDE_STAR_SELECT);
+                }
+
                 operationStack.push(GUIDE_CAPTURE);
 
                 // If we are being ask to go full frame, let's do that first
@@ -2029,7 +2036,18 @@ bool Guide::executeOperationStack()
 
         case GUIDE_CALIBRATING:
             if (guiderType == GUIDE_INTERNAL)
+            {
                 guider->setStarPosition(starCenter);
+                dynamic_cast<InternalGuider*>(guider)->setImageGuideEnabled(Options::imageGuidingEnabled());
+
+                // No need to calibrate
+                if (Options::imageGuidingEnabled())
+                {
+                    setStatus(GUIDE_CALIBRATION_SUCESS);
+                    break;
+                }
+            }
+
             if (guider->calibrate())
             {
                 if (guiderType == GUIDE_INTERNAL)

kstars/ekos/guide/internalguide/gmath.cpp

@@ -20,6 +20,7 @@
 #include "matr.h"
 
 #include "fitsviewer/fitsview.h"
+#include "imageautoguiding.h"
 
 #define DEF_SQR_0	(8-0)
 #define DEF_SQR_1	(16-0)
@@ -106,7 +107,10 @@ cgmath::~cgmath()
     delete [] drift[GUIDE_RA];
     delete [] drift[GUIDE_DEC];
 
+    foreach(float *region, referenceRegions)
+        delete [] region;
 
+    referenceRegions.clear();
 }
 
 
@@ -505,6 +509,20 @@ void cgmath::start( void )
     delta_prev = sigma_prev = sigma = 0;
 
     preview_mode = false;
+
+    // Create reference Image
+    if (imageGuideEnabled)
+    {
+        foreach(float *region, referenceRegions)
+            delete [] region;
+
+        referenceRegions.clear();
+
+        referenceRegions = partitionImage();
+
+        reticle_pos = Vector( 0, 0, 0 );
+    }
+
 }
 
 
@@ -535,6 +553,149 @@ void cgmath::setLostStar(bool is_lost)
     lost_star = is_lost;
 }
 
+float * cgmath::createFloatImage() const
+{
+    FITSData * imageData = guideView->getImageData();
+
+    // #1 Convert to float array
+    // We only process 1st plane if it is a color image
+    uint32_t imgSize = imageData->getSize();
+    float *imgFloat = new float[imgSize];
+    if (imgFloat == NULL)
+    {
+        qCritical() << "Not enough memory for float image array!";
+        return NULL;
+    }
+
+    switch (imageData->getDataType())
+    {
+    case TBYTE:
+    {
+        uint8_t *buffer = imageData->getImageBuffer();
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TSHORT:
+    {
+        int16_t *buffer = reinterpret_cast<int16_t*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TUSHORT:
+    {
+        uint16_t *buffer = reinterpret_cast<uint16_t*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TLONG:
+    {
+        int32_t *buffer = reinterpret_cast<int32_t*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TULONG:
+    {
+        uint32_t *buffer = reinterpret_cast<uint32_t*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TFLOAT:
+    {
+        float *buffer = reinterpret_cast<float*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TLONGLONG:
+    {
+        int64_t *buffer = reinterpret_cast<int64_t*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+        break;
+
+    case TDOUBLE:
+    {
+        double *buffer = reinterpret_cast<double*>(imageData->getImageBuffer());
+        for (uint32_t i=0; i < imgSize; i++)
+            imgFloat[i] = buffer[i];
+    }
+    break;
+
+    default:
+        return NULL;
+    }
+
+    return imgFloat;
+}
+
+QVector<float*> cgmath::partitionImage() const
+{
+    QVector<float*> regions;
+
+    FITSData * imageData = guideView->getImageData();
+
+    float *imgFloat = createFloatImage();
+
+    if (imgFloat == NULL)
+        return regions;
+
+    const uint32_t width      = imageData->getWidth();
+    const uint32_t height     = imageData->getHeight();
+
+    uint8_t xRegions = floor(width/regionAxis);
+    uint8_t yRegions = floor(height/regionAxis);
+    // Find number of regions to divide the image
+    //uint8_t regions =  xRegions * yRegions;
+
+    float *regionPtr = imgFloat;
+
+    for (uint8_t i=0; i < yRegions; i++)
+    {
+        for (uint8_t j=0; j < xRegions; j++)
+        {
+            // Allocate space for one region
+            float *oneRegion = new float[regionAxis*regionAxis];
+            // Create points to region and current location of the source image in the desired region
+            float *oneRegionPtr = oneRegion, *imgFloatPtr = regionPtr + j * regionAxis;
+
+            // copy from image to region line by line
+            for (uint32_t line=0; line < regionAxis; line++)
+            {
+                memcpy(oneRegionPtr, imgFloatPtr, regionAxis);
+                oneRegionPtr += regionAxis;
+                imgFloatPtr  += width;
+            }
+
+            regions.append(oneRegion);
+        }
+
+        // Move regionPtr block by (width * regionAxis) elements
+        regionPtr += width * regionAxis;
+    }
+
+    // We're done with imgFloat
+    delete [] imgFloat;
+
+    return regions;
+}
+
+void cgmath::setRegionAxis(const uint32_t &value)
+{
+    regionAxis = value;
+}
+
 Vector cgmath::findLocalStarPosition( void ) const
 {
     if (useRapidGuide)
@@ -544,6 +705,67 @@ Vector cgmath::findLocalStarPosition( void ) const
 
     FITSData * imageData = guideView->getImageData();
 
+    if (imageGuideEnabled)
+    {
+        float xshift=0, yshift=0;
+
+        QVector<Vector> shifts;
+        float xsum=0, ysum=0;
+
+        QVector<float *> imageParition = partitionImage();
+
+        if (imageParition.isEmpty())
+        {
+            qWarning() << "Failed to partiion regions in image!";
+            return Vector(-1,-1,-1);
+        }
+
+        if (imageParition.count() != referenceRegions.count())
+        {
+            qWarning() << "Mismatch between reference regions #" << referenceRegions.count() << "and image parition regions #" << imageParition.count();
+            // Clear memory in case of mis-match
+            foreach(float *region, imageParition)
+            {
+                delete [] region;
+            }
+
+            return Vector(-1,-1,-1);
+        }
+
+        for (uint8_t i=0; i < imageParition.count(); i++)
+        {
+            ImageAutoGuiding::ImageAutoGuiding1(referenceRegions[i], imageParition[i], regionAxis, &xshift, &yshift);
+            Vector shift(xshift, yshift, -1);
+            if (Options::guideLogging())
+                qDebug() << "Guide: Region #" << i << ": X-Shift=" << xshift << "Y-Shift=" << yshift;
+
+            xsum += xshift;
+            ysum += yshift;
+            shifts.append(shift);
+        }
+
+        // Delete partitions
+        foreach(float *region, imageParition)
+        {
+            delete [] region;
+        }
+        imageParition.clear();
+
+        float average_x= xsum / referenceRegions.count();
+        float average_y= ysum / referenceRegions.count();
+
+        float median_x = shifts[referenceRegions.count()/2-1].x;
+        float median_y = shifts[referenceRegions.count()/2-1].y;
+
+        if (Options::guideLogging())
+        {
+            qDebug() << "Guide: Average : X-Shift=" << average_x << "Y-Shift=" << average_y;
+            qDebug() << "Guide: Median  : X-Shift=" << median_x << "Y-Shift=" << median_y;
+        }
+
+        return Vector(median_x, median_y, -1);
+    }
+
     switch (imageData->getDataType())
     {
         case TBYTE:
@@ -1203,6 +1425,16 @@ const char * cgmath::get_direction_string(GuideDirection dir)
 
 }
 
+bool cgmath::isImageGuideEnabled() const
+{
+    return imageGuideEnabled;
+}
+
+void cgmath::setImageGuideEnabled(bool value)
+{
+    imageGuideEnabled = value;
+}
+
 //---------------------------------------------------------------------------------------
 cproc_in_params::cproc_in_params()
 {

kstars/ekos/guide/internalguide/gmath.h

@@ -189,7 +189,12 @@ class cgmath : public QObject
         // Logging
         void setLogFile(QFile * file);
 
-    signals:
+        bool isImageGuideEnabled() const;
+        void setImageGuideEnabled(bool value);
+
+        void setRegionAxis(const uint32_t &value);
+
+signals:
         void newAxisDelta(double delta_ra, double delta_dec);
         void newStarPosition(QVector3D, bool);
 
@@ -240,9 +245,19 @@ class cgmath : public QObject
         const char * get_direction_string(GuideDirection dir);
 
         // rapid guide
-        bool useRapidGuide;
+        bool useRapidGuide=false;
         double rapidDX, rapidDY;
 
+        // Image Guide
+        bool imageGuideEnabled=false;
+        // Creates a new float image from the guideView image data. The returned image MUST be deleted later or memory will leak.
+        float * createFloatImage() const;
+        // Partition guideView image into NxN square regions each of size axis*axis. The returned vector contains pointers to
+        // the newly allocated square images. It MUST be deleted later by delete[] or memory will leak.
+        QVector<float*> partitionImage() const;
+        uint32_t regionAxis=64;
+        QVector<float*> referenceRegions;
+
         // dithering
         double ditherRate[2];
 

kstars/ekos/guide/internalguide/imageautoguiding.cpp

+/*  Image Guide Algorithm
+    Copyright (C) 2017 Bob Majewski <cygnus257@yahoo.com>
+
+    This application is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public
+    License as published by the Free Software Foundation; either
+    version 2 of the License, or (at your option) any later version.
+ */
+
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "imageautoguiding.h"
+
+#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
+#define NR_END 1
+#define FREE_ARG char*
+
+
+#define	TWOPI	6.28318530717959
+#define FFITMAX	0.05
+
+
+//void ImageAutoGuiding1(float *ref,float *im,int n,float *xshift,float *yshift);
+
+float ***f3tensorSP(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_f3tensorSP(float ***t, long nrl, long nrh, long ncl, long nch,
+    long ndl, long ndh);
+
+float **matrixSP(long nrl, long nrh, long ncl, long nch);
+
+void free_matrixSP(float **m, long nrl, long nrh, long ncl, long nch);
+void nrerrorNR(void);
+
+void fournNR(float data[], unsigned long  nn[], long ndim, long isign);
+
+void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2,
+           unsigned long nn3, long isign);
+
+void ShiftEST(float ***testimage,float ***refimage,int n,
+float *xshift,float *yshift,int k);
+
+namespace ImageAutoGuiding
+{
+
+void ImageAutoGuiding1(float *ref,float *im,int n,float *xshift,float *yshift)
+{
+	float ***RefImage,***TestImage;
+	int i,j,k;
+	float x,y;
+    
+    /* Allocate memory */
+	
+	RefImage = f3tensorSP(1,1,1,n,1,n);
+	TestImage = f3tensorSP(1,1,1,n,1,n);
+	
+	/* Load Data */
+	
+	k = 0;
+	
+	for(j=1;j<=n;++j)
+	{
+	
+		for(i=1;i<=n;++i)
+		{
+			RefImage[1][j][i] = ref[k];
+			TestImage[1][j][i] = im[k];
+		
+			k += 1;
+		}
+	
+	}
+	
+	/* Calculate Image Shifts  */
+	
+	ShiftEST(TestImage,RefImage,n,&x,&y,1);
+	
+	*xshift = x;
+	*yshift = y;
+    
+    /* free memory */
+	
+	free_f3tensorSP(RefImage,1,1,1,n,1,n);
+	free_f3tensorSP(TestImage,1,1,1,n,1,n);
+	
+	
+	
+}
+
+}
+
+// Calculates Image Shifts
+
+
+void ShiftEST(float ***testimage,float ***refimage,int n,
+float *xshift,float *yshift,int k)
+{
+int  	ix,iy,nh,nhplusone;
+double	deltax,deltay,fx2sum,fy2sum,phifxsum,phifysum,fxfysum;
+double fx,fy,ff,fn,re,im,testre,testim,rev,imv,phi;
+double power,dem,f2,f2limit;
+float **speq;
+
+
+
+	f2limit = FFITMAX*FFITMAX;
+	
+	speq = matrixSP(1,1,1,2*n);
+	
+	
+	nh = n/2;
+	nhplusone = nh + 1;
+	
+	
+	fn = ((float) n);
+	ff = 1.0/fn;
+
+	
+	/* FFT of Reference */
+	
+	for(ix=1;ix<=2*n;++ix)
+	{
+		speq[1][ix] = 0.0;
+	}
+	
+	rlft3NR(refimage,speq,1,n,n,1);
+	
+	/* FFT of Test Image */
+	
+	for(ix=1;ix<=2*n;++ix)
+	{
+		speq[1][ix] = 0.0;
+	}
+	
+	rlft3NR(testimage,speq,1,n,n,1);
+	
+	/* Solving for slopes  */
+		
+	fx2sum = 0.0;
+	fy2sum = 0.0;
+	
+	phifxsum = 0.0;
+	phifysum = 0.0;
+	
+	fxfysum = 0.0;
+	
+	
+	for(ix =1;ix <= n;++ix)
+	{
+	
+		
+		if(ix <= nhplusone)
+		{
+			fx = ff*((float) (ix-1));
+		}
+		else
+		{
+			fx = - ff*((float) (n + 1 - ix));
+		} 
+			
+			
+				for(iy = 1;iy <= nh;++iy)
+				{
+				
+					fy = ff*((float) (iy-1));
+					
+					f2 = fx*fx + fy*fy;
+					
+                    /* Limit to Low Spatial Frequencies */
+					
+					if(f2 < f2limit)
+					{
+						
+				
+						/* Real part reference image */
+					
+						re = refimage[k][ix][2*iy-1];
+					
+						/* Imaginary part reference image */
+					
+						im = refimage[k][ix][2*iy];
+						
+						power = re*re + im*im;
+						
+						
+						/* Real part test image */
+					
+						testre = testimage[k][ix][2*iy-1];
+					
+						/* Imaginary part image */
+					
+						testim = testimage[k][ix][2*iy];
+					
+						
+						rev = re*testre + im*testim;
+						imv = re*testim - im*testre;
+                        
+                        /* Find Phase */
+						
+						phi = atan2(imv,rev);
+						
+												
+						fx2sum += power*fx*fx;
+						fy2sum += power*fy*fy;
+						
+						phifxsum += power*fx*phi;
+						phifysum += power*fy*phi;
+						
+						fxfysum += power*fx*fy;
+						
+						
+					
+					
+					}
+					
+				}
+			}
+			
+			
+			
+			
+			/* calculate subpixel shift */
+	
+			dem = fx2sum*fy2sum - fxfysum*fxfysum;
+			
+	
+            deltax = (phifxsum*fy2sum - fxfysum*phifysum)/(dem*TWOPI);
+            deltay = (phifysum*fx2sum - fxfysum*phifxsum)/(dem*TWOPI);
+			
+        
+			
+			free_matrixSP(speq,1,1,1,2*n);
+			
+			/* You can change the shift mapping here */
+	
+			
+			*xshift = deltax;
+			*yshift = deltay;
+		
+
+			
+}
+
+
+// 2 and 3 Dimensional FFT Routine for Real Data
+// Numerical Recipes in C Second Edition
+// The Art of Scientific Computing
+// 1999
+
+void rlft3NR(float ***data, float **speq, unsigned long nn1, unsigned long nn2,
+           unsigned long nn3, long isign)
+{
+    void fournNR(float data[], unsigned long nn[], long ndim, long isign);
+    void nrerror();
+    unsigned long i1,i2,i3,j1,j2,j3,nn[4],ii3;
+    double theta,wpi,wpr,wtemp;
+    float c1,c2,h1r,h1i,h2r,h2i;
+    float wi,wr;
+    
+    if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3)
+        nrerrorNR();
+    c1=0.5;
+    c2 = -0.5*isign;
+    theta=isign*(6.28318530717959/nn3);
+    wtemp=sin(0.5*theta);
+    wpr = -2.0*wtemp*wtemp;