IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 29, 2008, 10:43:18 AM (18 years ago)
Author:
eugene
Message:

adding a threaded version of psImageSmoothMask (only supports F32 for now)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageConvolve.c

    r20460 r20462  
    77/// @author Eugene Magnier, IfA
    88///
    9 /// @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
    10 /// @date $Date: 2008-10-29 19:46:31 $
     9/// @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
     10/// @date $Date: 2008-10-29 20:43:18 $
    1111///
    1212/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    711711}
    712712
    713 
    714 
    715713psImage *psImageSmoothMask(psImage *output,
    716714                           const psImage *image,
    … …  
    751749
    752750    switch (image->type.type) {
    753         IMAGESMOOTH_MASK_CASE(F32);
    754         IMAGESMOOTH_MASK_CASE(F64);
     751      case PS_TYPE_F32: {
     752          psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
     753          psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_MASK); /* Mask for calculation image; BW */
     754   
     755          /** Smooth in X direction **/
     756          for (int j = 0; j < numRows; j++) {
     757              for (int i = 0; i < numCols; i++) {
     758                  int xMin = PS_MAX(i - size, 0);
     759                  int xMax = PS_MIN(i + size, xLast);
     760                  const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[j][xMin];
     761                  const psF32 *imageData = &image->data.F32[j][xMin];
     762                  int uMin = - PS_MIN(i, size); /* Minimum kernel index */
     763                  const psF32 *gaussData = &gauss[uMin];
     764                  double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
     765                  for (int x = xMin; x <= xMax; x++, maskData++, imageData++, gaussData++) {
     766                      if (*maskData & maskVal) {
     767                          continue;
     768                      }
     769                      // if ((i == 1000) && (j > 10) && (j < 15)) {
     770                      //          fprintf (stderr, "%d %d  %d  %f  %f   %f %f\n", i, j, x, *imageData, *gaussData, sumIG, sumG);
     771                      // }
     772                      sumIG += *imageData * *gaussData;
     773                      sumG += *gaussData;
     774                  }
     775                  if (sumG > minGauss) {
     776                      /* BW */
     777                      calculation->data.F32[i][j] = sumIG / sumG;
     778                      calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0;
     779                  } else {
     780                      /* BW */
     781                      calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0xFF;
     782                  }
     783              }
     784          }
     785   
     786          output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
     787   
     788          /** Smooth in Y direction  **/
     789          for (int i = 0; i < numCols; i++) {
     790              for (int j = 0; j < numRows; j++) {
     791                  int yMin = PS_MAX(j - size, 0);
     792                  int yMax = PS_MIN(j + size, yLast);
     793                  const psMaskType *maskData = &calcMask->data.PS_TYPE_MASK_DATA[i][yMin]; /* BW */
     794                  const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
     795                  int vMin = - PS_MIN(j, size); /* Minimum kernel index */
     796                  const psF32 *gaussData = &gauss[vMin];
     797                  double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
     798                  for (int y = yMin; y <= yMax; y++, maskData++, imageData++, gaussData++) {
     799                      if (*maskData) {
     800                          continue;
     801                      }
     802                      sumIG += *imageData * *gaussData;
     803                      sumG += *gaussData;
     804                  }
     805                  output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
     806              }
     807          }
     808   
     809          psFree(calculation);
     810          psFree(calcMask);
     811          psFree(gaussNorm);
     812   
     813          return output;
     814      }
     815        IMAGESMOOTH_MASK_CASE(F64);
     816      default:
     817        psFree(gaussNorm);
     818        char *typeStr;
     819        PS_TYPE_NAME(typeStr,image->type.type);
     820        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     821                _("Specified psImage type, %s, is not supported."),
     822                typeStr);
     823        return false;
     824    }
     825
     826    psAbort("Should never reach here.");
     827    return false;
     828}
     829
     830bool psImageSmoothMask_ScanRows_F32 (psImage *calculation,
     831                                     psImage *calcMask,
     832                                     const psImage *image,
     833                                     const psImage *mask,
     834                                     psMaskType maskVal,
     835                                     psVector *gaussNorm,
     836                                     float minGauss,
     837                                     int size,
     838                                     int rowStart,
     839                                     int rowStop)
     840{
     841    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
     842
     843    int numCols = image->numCols;
     844    int xLast = numCols - 1; // Last index
     845
     846    for (int j = rowStart; j < rowStop; j++) {
     847        for (int i = 0; i < numCols; i++) {
     848            int xMin = PS_MAX(i - size, 0);
     849            int xMax = PS_MIN(i + size, xLast);
     850            const psMaskType *maskData = &mask->data.PS_TYPE_MASK_DATA[j][xMin];
     851            const psF32 *imageData = &image->data.F32[j][xMin];
     852            int uMin = - PS_MIN(i, size); /* Minimum kernel index */
     853            const psF32 *gaussData = &gauss[uMin];
     854            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
     855            for (int x = xMin; x <= xMax; x++, maskData++, imageData++, gaussData++) {
     856                if (*maskData & maskVal) {
     857                    continue;
     858                }
     859                sumIG += *imageData * *gaussData;
     860                sumG += *gaussData;
     861            }
     862            if (sumG > minGauss) {
     863                /* BW */
     864                calculation->data.F32[i][j] = sumIG / sumG;
     865                calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0;
     866            } else {
     867                /* BW */
     868                calcMask->data.PS_TYPE_MASK_DATA[i][j] = 0xFF;
     869            }
     870        }
     871    }
     872    return true;
     873}
     874
     875bool psImageSmoothMask_ScanCols_F32 (psImage *output,
     876                                     psImage *calculation,
     877                                     psImage *calcMask,
     878                                     psMaskType maskVal,
     879                                     psVector *gaussNorm,
     880                                     float minGauss,
     881                                     int size,
     882                                     int colStart,
     883                                     int colStop)
     884{
     885    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
     886
     887    int numRows = output->numRows;
     888    int yLast = numRows - 1; // Last index
     889
     890    for (int i = colStart; i < colStop; i++) {
     891        for (int j = 0; j < numRows; j++) {
     892            int yMin = PS_MAX(j - size, 0);
     893            int yMax = PS_MIN(j + size, yLast);
     894            const psMaskType *maskData = &calcMask->data.PS_TYPE_MASK_DATA[i][yMin]; /* BW */
     895            const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
     896            int vMin = - PS_MIN(j, size); /* Minimum kernel index */
     897            const psF32 *gaussData = &gauss[vMin];
     898            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
     899            for (int y = yMin; y <= yMax; y++, maskData++, imageData++, gaussData++) {
     900                if (*maskData) {
     901                    continue;
     902                }
     903                sumIG += *imageData * *gaussData;
     904                sumG += *gaussData;
     905            }
     906            output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
     907        }
     908    }
     909    return true;
     910}
     911   
     912bool psImageSmoothMask_ScanRows_F32_Threaded(psThreadJob *job)
     913{
     914    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     915    psAssert (job->args, "programming error: job args not alloced?");
     916    psAssert (job->args->n == 10, "programming error: job Nargs mismatch");
     917
     918    psImage *calculation  = job->args->data[0]; // calculation image
     919    psImage *calcMask     = job->args->data[1]; // calculation mask
     920    const psImage *image  = job->args->data[2]; // input image
     921    const psImage *mask   = job->args->data[3]; // input mask
     922
     923    psMaskType maskVal    = PS_SCALAR_VALUE(job->args->data[4],U8);
     924    psVector *gaussNorm   = job->args->data[5]; // gauss kernel
     925    float minGauss        = PS_SCALAR_VALUE(job->args->data[6],F32);
     926    int size              = PS_SCALAR_VALUE(job->args->data[7],S32);
     927    int rowStart          = PS_SCALAR_VALUE(job->args->data[8],S32);
     928    int rowStop           = PS_SCALAR_VALUE(job->args->data[9],S32);
     929    return psImageSmoothMask_ScanRows_F32(calculation, calcMask, image, mask, maskVal,
     930                                          gaussNorm, minGauss, size,
     931                                          rowStart, rowStop);
     932}
     933
     934bool psImageSmoothMask_ScanCols_F32_Threaded(psThreadJob *job)
     935{
     936    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     937    psAssert (job->args, "programming error: job args not alloced?");
     938    psAssert (job->args->n == 9, "programming error: job Nargs mismatch");
     939
     940    psImage *output       = job->args->data[0]; // input image
     941    psImage *calculation  = job->args->data[1]; // calculation image
     942    psImage *calcMask     = job->args->data[2]; // calculation mask
     943    psMaskType maskVal    = PS_SCALAR_VALUE(job->args->data[3],U8);
     944
     945    psVector *gaussNorm   = job->args->data[4]; // gauss kernel
     946    float minGauss        = PS_SCALAR_VALUE(job->args->data[5],F32);
     947    int size              = PS_SCALAR_VALUE(job->args->data[6],S32);
     948    int colStart          = PS_SCALAR_VALUE(job->args->data[7],S32);
     949    int colStop           = PS_SCALAR_VALUE(job->args->data[8],S32);
     950    return psImageSmoothMask_ScanCols_F32(output, calculation, calcMask, maskVal,
     951                                          gaussNorm, minGauss, size,
     952                                          colStart, colStop);
     953}
     954
     955psImage *psImageSmoothMask_Threaded(psImage *output,
     956                                    const psImage *image,
     957                                    const psImage *mask,
     958                                    psMaskType maskVal,
     959                                    float sigma,
     960                                    float numSigma,
     961                                    float minGauss)
     962{
     963    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
     964    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
     965    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL);
     966    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
     967
     968    int numCols = image->numCols, numRows = image->numRows; // Size of image
     969
     970    int nThreads = psThreadPoolSize();
     971    int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols;
     972    int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows;
     973
     974    // relevant terms
     975    int size = sigma * numSigma + 0.5;  // Half-size of kernel
     976    int fullSize = 2*size + 1;          // Total number of pixels in convolution kernel
     977
     978    // Generate normalized gaussian
     979    psVector *gaussNorm = psVectorAlloc(fullSize + 1, PS_TYPE_F32); // Normalised Gaussian
     980    {
     981        double sum = 0.0;
     982        double norm = -0.5/(sigma*sigma);
     983        for (int i = 0, u = -size; i <= fullSize; i++, u++) {
     984            float value = expf(norm*PS_SQR(u));
     985            gaussNorm->data.F32[i] = value;
     986            sum += value;
     987        }
     988        sum = 1.0 / sum;
     989        for (int i = 0; i <= fullSize; i++) {
     990            gaussNorm->data.F32[i] *= sum;
     991        }
     992    }
     993
     994    switch (image->type.type) {
     995        // Smooth an image with masked pixels
     996        // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
     997        // [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW"
     998
     999      case PS_TYPE_F32: {
     1000          psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
     1001          psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_MASK); /* Mask for calculation image; BW */
     1002   
     1003          /** Smooth in X direction **/
     1004          for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) {
     1005              int rowStop = PS_MIN (rowStart + scanRows, numRows);
     1006
     1007              // allocate a job, construct the arguments for this job
     1008              psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
     1009              psArrayAdd(job->args, 1, calculation);
     1010              psArrayAdd(job->args, 1, calcMask);
     1011              psArrayAdd(job->args, 1, (psImage *) image); // cast away const
     1012              psArrayAdd(job->args, 1, (psImage *) mask); // cast away const
     1013              PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
     1014              psArrayAdd(job->args, 1, gaussNorm);
     1015              PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
     1016              PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
     1017              PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32);
     1018              PS_ARRAY_ADD_SCALAR(job->args, rowStop,  PS_TYPE_S32);
     1019              // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, gauss, minGauss, size, rowStart, rowStop);
     1020
     1021              // if threading is not active, we simply run the job and return
     1022              if (!psThreadJobAddPending(job)) {
     1023                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
     1024                  psFree(job);
     1025                  return false;
     1026              }
     1027              psFree(job);
     1028
     1029          }
     1030
     1031          // wait here for the threaded jobs to finish (NOP if threading is not active)
     1032          if (!psThreadPoolWait(true)) {
     1033              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
     1034              return false;
     1035          }
     1036
     1037          output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
     1038   
     1039          /** Smooth in Y direction  **/
     1040          for (int colStart = 0; colStart < numCols; colStart+=scanCols) {
     1041              int colStop = PS_MIN (colStart + scanCols, numCols);
     1042
     1043              // allocate a job, construct the arguments for this job
     1044              psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
     1045              psArrayAdd(job->args, 1, output);
     1046              psArrayAdd(job->args, 1, calculation);
     1047              psArrayAdd(job->args, 1, calcMask);
     1048              PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_U8);
     1049              psArrayAdd(job->args, 1, gaussNorm);
     1050              PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
     1051              PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
     1052              PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32);
     1053              PS_ARRAY_ADD_SCALAR(job->args, colStop,  PS_TYPE_S32);
     1054              // -> psImageSmoothMask_ScanCols_F32 (output, calculation, calcMask, maskVal, gauss, minGauss, size, colStart, colStop);
     1055
     1056              // if threading is not active, we simply run the job and return
     1057              if (!psThreadJobAddPending(job)) {
     1058                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
     1059                  psFree(job);
     1060                  return false;
     1061              }
     1062              psFree(job);
     1063          }
     1064
     1065          // wait here for the threaded jobs to finish (NOP if threading is not active)
     1066          if (!psThreadPoolWait(true)) {
     1067              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
     1068              return false;
     1069          }
     1070          psFree(calculation);
     1071          psFree(calcMask);
     1072          psFree(gaussNorm);
     1073   
     1074          return output;
     1075      }
    7551076      default:
    7561077        psFree(gaussNorm);
    … …  
    10471368    if (threaded) {
    10481369        int numThreads = psThreadPoolSize(); // Number of threads
    1049         float cols = (float)numCols / (float)numThreads; // Number of cols to do at once
    1050         for (int i = 0; i < numThreads; i++) {
    1051             int start = i * cols;       // Starting colunms
    1052             int stop = (i + 1) * cols;  // Stopping columns
     1370        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
     1371        for (int start = 0; start < numRows; start += deltaRows) {
     1372            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
    10531373
    10541374            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
    … …  
    10601380            PS_ARRAY_ADD_SCALAR(job->args, xMin, PS_TYPE_S32);
    10611381            PS_ARRAY_ADD_SCALAR(job->args, xMax, PS_TYPE_S32);
    1062             PS_ARRAY_ADD_SCALAR(job->args, 0x00, PS_TYPE_U8);
     1382            PS_ARRAY_ADD_SCALAR(job->args, 0x00, PS_TYPE_U8); // specify rows
    10631383            if (!psThreadJobAddPending(job)) {
    10641384                psFree(job);
    … …  
    10861406    if (threaded) {
    10871407        int numThreads = psThreadPoolSize(); // Number of threads
    1088         float cols = (float)numCols / (float)numThreads; // Number of columns to do at once
    1089         for (int i = 0; i < numThreads; i++) {
    1090             int start = i * cols;       // Starting column
    1091             int stop = (i + 1) * cols;  // Stopping column
     1408        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
     1409        for (int start = 0; start < numCols; start += deltaCols) {
     1410            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
    10921411
    10931412            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
     1413            psArrayAdd(job->args, 1, out);
    10941414            psArrayAdd(job->args, 1, conv);
    1095             psArrayAdd(job->args, 1, out);
    10961415            PS_ARRAY_ADD_SCALAR(job->args, start, PS_TYPE_S32);
    10971416            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
    … …  
    10991418            PS_ARRAY_ADD_SCALAR(job->args, yMin, PS_TYPE_S32);
    11001419            PS_ARRAY_ADD_SCALAR(job->args, yMax, PS_TYPE_S32);
    1101             PS_ARRAY_ADD_SCALAR(job->args, 0xff, PS_TYPE_U8);
     1420            PS_ARRAY_ADD_SCALAR(job->args, 0xff, PS_TYPE_U8); // specify cols
    11021421            if (!psThreadJobAddPending(job)) {
    11031422                psFree(job);
    … …  
    11861505}
    11871506
    1188 
     1507// XXX for now, either thread all or none
     1508// have to call this before calling psImageSmoothMask_Threaded
    11891509void psImageConvolveSetThreads(bool set)
    11901510{
    … …  
    11961516            psFree(task);
    11971517        }
     1518        {
     1519            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANROWS", 10);
     1520            task->function = &psImageSmoothMask_ScanRows_F32_Threaded;
     1521            psThreadTaskAdd(task);
     1522            psFree(task);
     1523        }
     1524        {
     1525            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS", 9);
     1526            task->function = &psImageSmoothMask_ScanCols_F32_Threaded;
     1527            psThreadTaskAdd(task);
     1528            psFree(task);
     1529        }
    11981530    } else if (!set && threaded) {
    11991531        psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");
    1200     }
    1201 
     1532        psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
     1533        psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
     1534    }
    12021535    threaded = set;
    1203 
    12041536}
    12051537
Note: See TracChangeset for help on using the changeset viewer.